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Preparation of Optical-Fiber Ends for 
Low-Loss Tape Splices 


By E. L. CHINNOCK, D. GLOGE, P. W. SMITH, and D. L. BISBEE 
(Manuscript received September 18, 1974) 


We describe a reliable method of preparing planar fiber tape ends by 
fiber fracture. Using this technique, with suitable precautions to preserve 
cleanliness during splice preparation, we have measured a splice loss of 
less than 0.25 dB in 99 percent of all attempts. 


The alignment of fibers in prefabricated grooves! so far remains the 
simplest and most reliable method of connecting fibers, even in the 
case of fiber cable subgroups (tapes).?? Two problems were identified as 
most serious: 


(t) The preparation of satisfactory fiber ends was found difficult 
in the case of tapes and cables, because all ends must be in one 
cross-sectional plane. Grinding and polishing has proven feasible, 
but may not be entirely satisfactory, particularly in the prep- 
aration of field splices. 

(it) The splice losses have been higher than expected on the basis 
of single-fiber splice tests and have been scattered over a wide 
range. 


Although we are not certain that groove alignment necessarily 
provides the best splicing technique, we have used an advanced form 
of this technique developed by Cherin’ to take a closer look at the 
problems identified above. Even if other techniques prove more 
promising later on, the problems mentioned may still be present in 
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some form or other and seem serious enough to require thorough 
analysis now. 

The preparation of fiber ends discussed here is a modification of 
fracture techniques reported earlier for single-fiber splices.4 The de- 
vice used for this purpose was a compact and simple hand tool that 
could easily be operated in a cramped and narrow space. The essential 
element of this tool was a spring-steel strip over which the fibers were 
stretched. The mechanical characteristics of this strip primarily 
determined the stress distribution in the fibers, and, thus, by a proper 
choice of strip thickness, we were able to choose the appropriate ratio 
of bending to tensile stress for the particular fibers to be fractured. 

To prepare the ends of the fiber tape for splicing, one proceeds as 
follows: 


(t) The plastic of the tape is removed over a short distance so 
that the fibers are exposed in the area where the end is to be 
prepared. 

(it) The tape is placed between a spring-steel strip and two fric- 
tion plates, so that the exposed area is located under a diamond 
stylus. 

(zit) The spring-steel strip and the tape are bent. At the same time 
the friction plates slide a small distance along the spring-steel 
strip. This sliding action exerts an additional amount of tension 
on the tape, so that the optimal ratio between longitudinal 
and bending stress is obtained in the fibers. 

(iv) The diamond stylus (tip radius 50 um) is now drawn across 
the exposed fibers to produce scores. The slight pressure of a few 
grams imparted by a phosphor bronze spring suffices to pro- 
duce scores of a few micrometers in depth. As each fiber is 
scored, a fracture starts at the score and proceeds across the 
fiber producing a flat surface perpendicular to the fiber axis. 


For a more detailed explanation of this process, see Ref. 4. The order 
of the steps explained above is not imperative. As an alternative, the 
scoring can be done before tension is applied. 

Figure 1 shows a typical array of fiber ends obtained in this way. 
The fibers used were multimode fibers having a high-silica core with a 
diameter of about 80 wm and an outer diameter of 120 um. To make a 
tape splice, we prepared the ends of two tapes in the way discussed 
above. The tapes were then positioned slightly above a small grooved 
chip made from lead, copper, or aluminum (see Fig. 2). The chip was 
roughly 1 cm long and 3 to 4 mm wide. The grooves were embossed 
using a stainless-steel head that had six adjacent 90-degree grooves, 
each 80 wm deep. As shown in Fig. 2, the fibers were lowered into 
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Fig. 1—Tape end prepared by simultaneous fracture of fibers. Fibers were epoxied 
together after end preparation to keep them aligned for electron micrograph process. 
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Fig. 2—Sketch of splice arrangement before assembly. 
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grooves by pressing from the top with a stack made of (from top to 
bottom) a lead weight, a foam rubber pad, and a sheet of gum rubber 
about 150 wm thick. After the fibers were lowered into the grooves, the 
tapes were gently pushed together in axial direction. 

Tapes 2 m long were used to measure the splice loss. To simulate 
longer lengths, we injected light with a power distribution approximat- 
ing the steady-state distribution of these particular fibers. We first 
made a large number of loss measurements in unbroken and unspliced 
tapes, switching back and forth between the five fibers of each tape 
to determine the measuring uncertainty. We found a distribution that 
had an rms value of 1.7 percent. The splice loss was then determined 
by measuring the transmission before and after a small part (a few 
centimeters) was removed from the middle of each tape and the ends 
spliced together as explained above, adding a drop of index-matching 
oil or glycerin before covering the arrangement with the gum rubber 
sheet. The optical loss of the length of fiber removed was insignificant. 
Figure 3 shows a histogram of the splice losses measured in 60 attempts. 
Evidently some loss values were negative as a result of the measuring 
inaccuracy. Figure 4 shows the (smoothed) cumulative loss distribu- 
tion as measured and after the 1.7-percent rms measuring uncertainty 
was discounted. Of all measurements, 99 percent show a loss of less 
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Fig. 3—Histogram of measured splice loss. 
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Fig. 4—Cumulative distribution of splice loss as measured and after discounting 
the 1.7-percent rms measuring uncertainty. 


than 0.25 dB. However, a large part of the scatter of these measure- 
ments results from the measuring process and the actual cumulative 
distribution would predict that 99 percent of all splices would have a 
loss of less than 0.15 dB. 

We attribute the low-loss values obtained not only to the quality of 
the end faces, but also to the extreme care taken during these measure- 
ments to keep the splice area clean. Earlier microscopic observations 
of groove-aligned fiber array splices taught us that, generally, losses 
in the 10-percent range can be correlated with a contamination in the 
splice area which reveals itself by large-angle scattering at the joint. 
In most cases, we were able to identify the contaminating material 
on the end face of the fiber even after the splice was taken apart ; these 
materials tenaciously adhere to the fiber surface often even after 
ordinary cleaning procedures. We learned that an extended period of 
ultrasonic cleaning with isopropyl] alcohol of all parts involved in the 
splice was necessary before the splice loss decreased to the levels 
measured. The contaminant is usually not added during the end prep- 
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aration; it is not the result of dust accumulation from the surrounding 
air caused, for example, by electrostatic forces. We believe that the 
contamination results from a contact of the fiber ends with con- 
taminated surfaces, such as the grooved chip or the rubber sheet. We 
believe, also, that this sensitivity to contamination is a sufficient reason 
to consider splicing processes in which the fiber end surfaces are 
prepared and exposed after alignment has been achieved.® 

Permanent splices on the basis of the techniques described here 
were prepared by replacing the index-matching oil with a special 
epoxy. This epoxy flows down along the grooves and the fibers and 
permanently attaches the tape ends to the embossed metal chip and 
the gum rubber sheet. The rubber sheet ends extending beyond the 
chip (see Fig. 2) are then folded around the chip and attached to its 
bottom surface. Figure 5 shows a finished five-fiber splice. Splices of 
this kind were found to have sufficient intrinsic strength to be used as 
splices of cable subgroups; additional armor would of course in this 
case be provided around a stack of such subgroup splices in a cable. 
The loss distribution for these permanent splices showed no deviation 
from that of splices made using index-matching fluid, and no aging 
effect was noticed, at least not within the period of a few days. The 
epoxy has a room-temperature curing time of several hours, but faster- 
curing epoxies are being studied and should replace the one used with- 
out significant alteration of the results. 





Fig. 5—Finished five-fiber splice. 
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We are grateful to R. D. Standley who prepared the electron micro- 
graph of Fig. 1. 
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All-Glass Optical-Fiber Tapes 


By D. L. BISBEE and P. W. SMITH 
(Manuscript received September 18, 1974) 


We propose and demonstrate a new approach to the problem of splicing 
optical fibers in a fiber cable. The optical fiber cable subgroups (tapes) 
are made in such a way that the relative positions of the optical fibers are 
accurately maintained. By using glass as a rigid matrix material in 
which the optical fibers are held, we demonstrate that a simple scoring and 
stressing technique can be used to simultaneously prepare all the fiber 
ends for splicing. 


The potential of optical fibers as transmission media for optical 
communications systems has stimulated much work on the various 
problems that need to be overcome before a practical system can be 
built. One of these problems involves the development of techniques 
for connecting and splicing these fibers and fiber cables. Although 
several laboratory techniques for splicing fibers and groups of fibers 
with low splice losses have now been developed, they are all relatively 
complex techniques that require operations of high precision and are 
thus difficult to carry out in the field. 

In this paper, we propose a different approach to the splicing 
problem. Linear arrays (‘‘tapes’’) of optical fibers have been suggested 
as building blocks for optical fiber cables,! and a number of techniques 
for producing plastic-bonded tapes have been investigated.? We propose 
here a technique for fabricating a precision all-glass optical fiber tape 
that would have considerable advantages with regard to splicing opera- 
tions. The precision operations would be performed during the manu- 
facturing of the tapes, and splicing operations in the field would become 
relatively simple. The basic idea is to make a fiber tape in which the 
optical fibers are held in a rigid matrix with their relative positions 
accurately maintained. Further, by using a glass for the rigid matrix 
material, we can greatly simplify the problem of preparing the optical 
fiber ends for splicing: All the fibers comprising the tape can be pre- 
pared for splicing in a single operation by utilizing the scoring and 
stressing technique described earlier in Ref. 3. 
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We made all-glass tapes by fusing conventional clad soda-lime- 
silicate glass optical fibers together with lower melting point glass 
fibers in a precision jig. To get a stable bond, the glasses must have 
very nearly the same thermal expansion coefficient even though their 
softening temperatures are different. Glasses with these characteristics 
can be obtained. 

As an example, let us consider glass composed of SiOz, Na.O, and 
CaO (“‘soda-lime-silicate”’ glass). Morey‘ defines the softening point 
of glass as that point at which the viscosity becomes 107-* poises. From 
Ref. 4, we see that by changing the composition of our soda-lime- 
silicate glass, it is possible to make one composition that has a vis- 
cosity of 107-° poises at 650°C and another with a viscosity two orders 
of magnitude greater than this at the same temperature. The thermal 
expansion coefficient of the first will be 9.7 X 10~® per °C, while that 
of the second will be 11.3 X 10-® per °C. Thus, we can obtain glasses 
that have viscosities different by two orders of magnitude at the 
softening point, while their thermal expansion coefficients differ by 
only 15 percent. By using such glasses, we could ensure that the de- 
formation during the fusing operation would take place only in the 
low-melting-point glass, and the optical fibers would be essentially 
undistorted. 

The optical fibers used for these experiments were conventional clad 
multimode fibers made from soda-lime-silicate glass. The low-melting- 
point fibers were made from ferrite sealing glass obtained from the 
Corning Glass Works (No. 8463). The optical fibers had a softening 
temperature of 650°C and a thermal expansion coefficient of 
9.7 X 10-* per °C, while the softening temperature of the low-melting- 
point fibers was = 375°C, and the thermal expansion coefficient was 
10.4 X 10-° per °C. 

Figure 1 is a schematic view of the precision guide used to fabricate 
the all-glass fiber tape. The optical fibers are held accurately in place 
by the precision guide, while the low-melting-point triangular fibers are 
introduced in such a way that they press against two adjacent optical 
fibers. By means of a heating element, the glass is fused at the contact 
points and then allowed to cool as the completed tape is pulled out of 
the guide. 

The heating element is a 500-um outer diameter nichrome wire 
mounted on a manipulator. When the heating wire is brought to 
within 100 um of the cold fibers and heated by passing current through 
it, the triangular fibers melt and form beads of molten glass that touch 
the heater. This glass is then smoothly spread along the round fibers 
as the fibers are fed into the guide, giving a uniform bond. 

The precision guide is made of boron nitride—a material that 
resembles soapstone. It is used because it is easily machinable to 
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Fig. 1—Schematic diagram of the device used to fuse the all-glass fiber tapes. 


close tolerances, has a slick surface, and is not affected by the heat 
applied. 

To feed the fibers through the precision guide, rollers are used as 
shown in Fig. 1. One pair of rollers is used to push the triangular fibers 
into contact with the round fibers, and one pair is used to pull the 
fused tape out of the guide. The lower roller in each pair is a 6.25-mm 
shaft with a tight-fitting sleeve of vinyl 0.5-mm thick giving an outer 
diameter of 7.25 mm. The upper roller of each pair is a 6.25-mm shaft 
with a flexible plastic sleeve giving an outer diameter of 9.38 mm. One 
of the smaller rollers is driven by a variable speed motor and is con- 
nected to the other small roller through an idler gear. Before entering 
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the guide, the fibers were cleaned by passing through a solvent-soaked 
wick. 

To make the tape, the heater current is turned on while the heating 
wire is still far from the fibers (more than 250 yum), the driving motor 
is turned on, and then, while the fibers are passing through the guide, 
the heating wire is lowered by a manipulator (while the operator 
watches through a microscope) until the triangular fibers melt and 
bond smoothly to the round fibers. The tape was formed at about 3 em 
per minute. Figure 2 is a photograph of the resultant fiber tape. 

It has previously been shown that, by scoring an optical fiber and 
subjecting it to a properly tailored stress distribution, a smooth 
fracture perpendicular to the fiber axis can be obtained.’ It has also 
been observed that such fracture behavior can be obtained with more 
complex fiber cross sections. We used a diamond stylus to score the 
fiber tape and fractured it by subjecting it to bending and tension 
stress. The fracturing operation was performed with the aid of a 
device similar to that described in Ref. 5. Internal strains would some- 
times cause the tape to fracture in such a way that the break was not 
perpendicular to the tape axis. Nevertheless, good fiber ends were 
usually produced. Figure 3 shows a fiber tape end produced in this way. 

The solder-glass fibers that were used tended to break very easily 





Fig. 2—A section of all-glass fiber tape. 
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Fig. 3—A fiber tape end prepared for splicing by the scoring and stressing technique. 


and had to be handled with great care; thus, the tapes we made were 
fragile. When fused together, the round and triangular fibers are so 
intimately bonded together that a crack originating in the glass of a 
triangular fiber would propagate across the whole tape. If low-melting- 
point glass with better mechanical properties were used, we would 
expect the finished tape to be appreciably stronger. As can be seen in 
Fig. 3, the tape is no thicker than its round fibers, so it is flexible in the 
direction of its thinner dimension. 

Because the optical fibers are accurately positioned during the 
manufacture of the tape, the most difficult part of the splicing problem 
has already been solved, and the splicing of these tapes merely in- 
volves preparing the tape ends by scoring and bending and placing 
the prepared ends in a suitable holder with matching fluid and a cover 
to hold the tape ends in place, as shown in Fig. 4. To make a permanent 
splice, a transparent index-matching epoxy may be used. Note that at 
no stage during the entire splicing process does one have to deal with 
single optical fibers, and that the tapes can be placed by hand in the 
splicing holder without the need for any precise visual alignment. 
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Fig. 4—Schematic representation of tape splicing operations. 


We would like to thank A. R. Tynes for pulling the optical fibers 
and the triangular solder glass fibers used for these experiments. 
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Synthesis of Speech From a Dynamic Model 
of the Vocal Cords and Vocal Tract 


By J. L. FLANAGAN, K. ISHIZAKA, and K. L. SHIPLEY 
(Manuscript received July 17, 1974) 


We describe a computer model of the human vocal cords and vocal tract 
that is amenable to dynamic control by parameters directly identified in 
the human physiology. The control format consequently provides an 
efficient, parsimonious description of speech information. The control 
parameters represent subglottal lung pressure, vocal-cord tension and rest 
opening, vocal-tract shape, and nasal coupling. Using these inputs, we 
synthesize vowel-consonant-vowel syllables to demonstrate the dynamic 
behavior of the cord/tract model. We show that inherent properties of the 
model duplicate phenomena observed in human speech; in particular, 
cord/tract acoustic interaction, cord vibration, and tract-wall radiation 
during occlusion, and voicing onset-offset behavior. Finally, we describe 
an approach to deriving the physiological controls automatically from 
printed text, and we present sentence-length synthesis obtained from a 
preluminary system. 


l. INTRODUCTION 


Speech sounds can be synthesized by a variety of means used to 
construct signal waveforms. Many ingenious methods have been re- 
corded. But speech synthesis generally has the practical purpose of 
producing intelligible sounds from control data that are as parsimonious 
as possible. In other words, the control data should represent an 
efficient, concise coding of the speech information. This motivation 
applies as much to analysis/synthesis techniques for speech transmis- 
sion as to computer voice-response systems which strive for efficient 
vocabulary storage and high versatility in message fabrication. 

Because speech is a human-generated signal, it is unlikely that a 
synthesis method can achieve the ultimate parsimony of input control 
without considerable attention to the parameters a human overtly 
manipulates in speaking. That is, one increases the information 
“built into” the synthesizer when its design exploits fundamental 
properties of the human speech mechanism. 

We therefore have chosen an approach to synthesis with which we 
can identify overtly the significant physiological parameters important 
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in speech production. Major system components obviously are the 
mechanism of voiced-sound generation and the mechanism for intel- 
ligibly modulating sound timbre, that is, the vocal cords and the vocal 
tract. Our approach, unlike that found conventionally in the speech 
literature, is not to make a linear separation of the sound source and 
vocal tract. More than this, we believe that source/tract interaction 
actually contributes built-in natural behavior that is significant in 
synthesis. This natural interaction is missing in approaches that 
assume linear separation of source and tract (unless provided at addi- 
tional expense and coding effort in the input data). 

The initial results stemming from this approach to synthesis are 
described below. 


Il. ACOUSTIC MODEL OF VOCAL CORDS AND VOCAL TRACT 


We view the acoustic system of the human vocal cords and vocal 
tract as shown at the top of Fig. 1. The lungs are an air reservoir, 
maintained at subglottal air pressure P, by contraction of the rib-cage 
muscles. The subglottal pressure is applied via the bronchi and trachea 
passages to the variable-area orifice controlled by the vocal cords. 

We model the cords as an acoustic-mechanical oscillator, wherein a 
single vocal cord is described by two masses, each having an associated 
stiffness and loss, which are “internally” coupled by a third stiffness. 
In previous work,!~* we established the philosophy leading to this 
description and gave a quantitative analysis of the vocal cord model. 

Oscillation of the vocal cord model results in the glottal volume 
velocity U,. This quantity typically has an impulsive waveform and it 
is the excitatory source for voiced sounds. 

The vocal tract proper is a nonuniform tube, about 17 cm long in 
man, extending from the cords to the mouth. Its cross-sectional area 
varies from zero to upwards of 20 cm?. The nasal tract is an ancillary 
tube about 60 cm! in total volume and coupled to the vocal tract by 
the trap-door action of the velum. Sound is radiated from the system 
as a result of the volume velocities at the mouth Um and nostril Un, 
and from vibration of the yielding sidewalls of the vocal tract. 

Cross-dimensions of the acoustic system are small compared to sound 
wavelengths of interest, and hence we confine our analysis to plane- 
wave propagation in the tract. We therefore represent the acoustic 
system as the bilateral, time-varying transmission line shown in the 
lower part of Fig. 1. Formulation of this system follows that given by 
Flanagan.® 

As illustrated, the lossy lung volume is ‘‘charged”’ to subglottal lung 
pressure P,, which is applied via the trachea-bronchi network, to the 
glottal (vocal-cord opening) impedance Z,. This nonlinear glottal 
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Fig. 1—Schematic diagram of the vocal cord/vocal tract system. 


impedance depends upon the glottal flow and area A,, which in turn 
depend upon the self-oscillating properties of the vocal cord model 
described in detail in an earlier paper.! The resulting volume flow U, 
is the excitation source for the vocal and nasal tracts. 

The shape of the vocal tract is defined by its cross-sectional area 
as a function of distance A(x), and the coupling to the time-invariant 
nasal cavity is governed by the velar impedance Z,. Volume velocity 
at mouth U,, and nostril U, flow through their respective radiation 
impedances Z, and Z,, both of which are in series with batteries 
representing the constant atmospheric pressure P,. (This formulation 
permits simulation of respiration as well.) The mouth and nostril 
radiation impedances are those for a circular piston in an infinite 
baffle.® 

Parameters of control for the speech synthesis system are the 
physiologically-based functions shown in Fig. 1. All vary with time. 
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They are subglottal lung pressure P,, vocal-cord tension Q, rest (or 
neutral) area of cord opening A,., nasal coupling N, and cross-sectional 
area function of the tract shape A(x). We are concerned here only with 
nonnasal sounds, hence nasal coupling will not figure in the discussion. 

Each T-section of the vocal-tract transmission line is represented in 
Fig. 2.5 An elemental length Az of the vocal tube has cross-sectional 
area A, terminal sound pressures p, and p2, and terminal volume veloci- 
ties U; and U». The sidewall has noninfinite mechanical impedance, 
and vibrates in response to the enclosed sound pressure with displace- 
ment £ This displacement radiates a per-unit-length sound pressure 
Pwa. Relations between terminal values of pressure and volume 





——» Up 





a L 
Pout = [P MOUTH *P NOSE le PWALL dx] 


(b) NETWORK 


Fig. 2—Circuit representation of plane acoustic wave propagation in an elemental 
length of tube with yielding sidewalls. 
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velocity for plane-wave propagation are given in the circuit in Fig. 2b, 
in which L is the per-unit-length inertance of the air mass of the tube 
element, & the viscous loss at the sidewall, G the heat-conduction loss, 
C the acoustic compliance of the contained air volume, Z, the acoustic 
equivalent mechanical impedance of the yielding wall, and Z, the 
radiation impedance of the wall, assumed to be that for a pulsating 
right circular cylinder.® 

The total sound output from the model is, following the long-wave 
assumptions, the linear superposition of the mouth and nostril radia- 
tion plus the spatially summated wall radiation. 

In addition, every T-section of the transmission-line network in- 
cludes a means for introducing turbulent noise excitation. This 
capability is provided by a series random pressure source Py with its 
internal resistance Ry, as shown in Fig. 3. This technique has been 
given in detail previously.4 The intensity (or rather mean-square 
variance) of the random pressure source is controlled by the Reynolds 
number of the flow at each network section, while the internal resis- 
tance is similarly modulated according to the Bernoulli loss in a 
constriction.® In both instances, the specified value of cross-sectional 
area A and the calculated resulting volume velocity flowing through the 
serial source completely describe the control functions. That is, no 
additional input data are required. 

More specifically, to simulate the conditions of turbulent-source 
generation in any section, the amplitude of the noise pressure is made 
directly proportional to the squared Reynolds number in excess of a 


L R R L r7x Rn 
ae ee ee ae RO 
Lw 
Rw 
G Cc 
Cw 
ZRw 


Fig. 3—Circuit representation of the turbulent noise source for each network 
section. The intensity of the random pressure source, Py, and its self-resistance, Ry, 
are controlled by the volume velocity and cross-sectional area at each section. 
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critical (threshold) Reynolds number for turbulent flow.4 The squared 
Reynolds number is proportional to U?/A, whereas the internal resist- 
ance of the turbulent source is, to first order, a flow-dependent loss 
proportional to |U|/A?. Therefore, as the prescribed section area 
becomes small in the presence of a large flow-velocity, turbulence condi- 
tions are indicated and the intensity of the noise source and the value 
of its internal resistance are increased. In the simulation, values of 
every dependent variable are calculated on a sample-by-sample basis 
to construct the time functions for the output sound pressure and all 
other pressure and velocity quantities.’ 

As a consequence of continually noting the magnitude of the cal- 
culated volume flow in each section and having the tract cross-sectional 
areas continually prescribed as input data, the synthesizer auto- 
matically introduces random noise excitation in any section when the 
Reynolds number is sufficiently high to indicate turbulent flow. The 
synthesizer, therefore, requires no additional data to produce voiceless 
sounds, but uses exactly the same control parameters to generate both 
voiced and voiceless sounds (or combinations of voiced and voiceless 
sounds). As Fig. 1 has shown, these control parameters are P,, Q, 
Ago, N, and A(z). 

As a practical matter in the computer implementation, we use a 
Py source produced from gaussian noise (or, rather, gaussian numbers) 
bandpass-filtered from 500 to 4,000 Hz. Further, to insure stability, 
the volume flow which modulates the serial noise source is low-pass 
filtered to 500 Hz. In other words, the noise source is modulated by 
low-frequency components of U, including the de flow. 

The transmission line model of Fig. 1 is described by a set of linear 
and nonlinear differential equations in which all coefficients also vary 
with time. This set of differential equations is approximated by 
difference equations, as previously described,? and programmed in a 
laboratory computer for on-line control. Twenty network sections are 
used to approximate the vocal tract. This formulation has permitted 
initial experiments with physiologically-based control of the synthesis 
model. 


Ill. ASSESSMENT OF WALL IMPEDANCE AND EFFECTS ON FORMANT 
BANDWIDTH 
All elements of the transmission line network have been well estab- 
lished in previous work, with the exception of the wall-vibration shunt 
branch of the circuit in Fig. 2. 
Assessment of wall effects in earlier calculations® utilized the only 
available mechanical impedance measurements of human tissue, 
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namely, chest, stomach, and thigh tissue. These data led to correct 
order-of-magnitude values for wall-vibration damping of formant 
resonances, but the values were clearly on the high side. 

To better assess the wall impedance, we have done two things. First, 
we have used data on human formant bandwidths to estimate contribu- 
tions to losses in our model. And second, we have made direct measure- 
ment of the mechanical impedance of the vocal-tract wall.” 

Formant bandwidths have been measured for the human vocal 
tract by van den Berg,*® Bogert,® Fujimura and Lindquist,!° and Dunn." 
Our programmed model allows us to adjust values of the wall-im- 
pedance parameters to effect three consistencies. It permits us to 
(t) adjust the wall-loss component to match glottis-closed formant 
bandwidths, (77) adjust the inductive reactance of the wall to produce 
the observed mouth-closed, lowest value of first formant frequency of 
about 200 Hz, and (277) choose a wall compliance to produce wall 
resonance substantially below 100 Hz. Small-signal-driven vibration of 
the cord oscillator in the model permits calculations of model response 
at any prescribed frequency. Furthermore, formant bandwidths mea- 
sured on real speech allow additional cross-checks of parameters used 
in the model formulation, especially for the loss components of the 
cord-oscillator source. Application of this knowledge in our model 
yields the formant bandwidth behavior shown in Fig. 4. 

In particular, Fig. 4 illustrates how the wall viscous loss parameter 
can be chosen to match glottis-closed formant bandwidth. This tech- 
nique has recently been analyzed in extensive quantitative form by 
Sondhi.” The wall loss is selected to match measured formant band- 
widths at formant frequencies around 300 to 500 Hz. In this fre- 
quency range, the contributions to formant bandwidth are mainly 
wall loss and glottal source loss. Viscosity, heat conduction at the 
walls, and mouth radiation resistance represent relatively small 
values (see Ref. 5, for example, for these calculation techniques). 
Note, too, that in Fig. 4 the vertically-sloping line of calculated band- 
width indicates the effect of wall impedance on the tuning of formant 
frequency. In the absence of additional data, we assume a uniform 
distribution of the per-unit-area wall impedance along the tract. The 
value we use for the mechanical per-unit-area impedance is 
(1600 + 71.5w)g/s/em?, where w is the radian frequency. This value is 
confirmed well by our direct measurements of wall impedance.’ 

Formant bandwidths measured in real speech"! permit a cross-check 
of the glottal oscillator parameters chosen in previous work.! Figure 4 
shows that the contribution to formant damping of the glottal source 
falls into the correct range of real speech measurements. This is a 
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_Fig. 4—Variation of formant bandwidth (damping) with formant frequency. The 
diagram shows the relative contributions of loss to the total formant bandwidth 
produced by the cord/tract model. 


gratifying confirmation of cord parameters chosen strictly on other 
bases, namely, according to physiological properties and oscillatory 
behavior.! 

The loss contributions of the glottal source in Fig. 4 are calculated 
for a nominal, midrange value of glottal rest area, namely A,. = 0.05 
em. The glottal contribution to formant bandwidth is, of course, a 
function of A,o. Figure 5 shows glottal loss contributions for other 
values of A, .. Note, especially, how the articulatory configuration 
of the tract influences the contribution of the glottal source to formant 
damping. 


IV. DYNAMIC BEHAVIOR OF THE CORD/TRACT MODEL 


How does this physiologically-based model of the vocal cords and 
vocal tract behave under dynamic control? Time-varying control 
inputs in the present study are P,, Q, A, and A(x). An obvious 
major problem is the determination of realistic values of these param- 
eters. As a first cut, fairly realistic data can be derived from direct 
measurements of lung pressure during speech," laryngeal muscle 
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Fig. 5—Variation of the glottal loss contribution to formant bandwidth as a 
function of formant frequency. The parameter is the neutral glottal area, Ago. 


electromyography,'* glottal transillumination,“ glottal pulses, and 
cine X-rays of the vocal tract.!® An important element at present is that 
all these data are not simultaneously available for a given subject. 
Experiments now under way aim to provide some simultaneous 
measurements.!” 

The physiological literature does provide adequate bases for dynamic 
tests on some simple utterances, using idealized input controls. We 
have, therefore, made first tests on vowel-consonant-vowel syllables 
(v-c-v) in which stress may be on either initial or final vowel, and 
where the intervocalic consonant is a voiced or unvoiced labial stop. 
These combinations also provide a convenient vehicle for exposing 
other physiologically realistic properties of the cord/tract model. 

Figure 6 shows the synthesis of the syllable /’aba/. Input controls 
are indicated in the top three traces. Because of the initial stress, Ps 
falls during the labial closure to a lower value. Because the intervocalic 
stop is voiced, A,, is maintained in a position favorable to cord oscilla- 
tion throughout. Cord tension, not shown, is also maintained constant. 
Any pitch changes are effected solely by P, variation and by the inter- 
action of tract load on the cord oscillator. Articulatory shape A (2) 
changes from /a — b — 2/. Because of space limitation in the illustra- 
tion, only the mouth area, Am, is displayed. 

Response behavior of the model to these input controls is shown in 
the bottom five traces: the sound spectrogram of the total output 
sound; A,; U,; the pressure waveform of the total output sound P; 
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Fig. 6—Control functions and sound output from the cord/tract model synthesizing 
the syllable /’aba/. The effects of sound radiation from the yielding sidewalls is 
evidenced in synthesized sound, and the vibration of the mouth cavity wall is illus- 
trated by the AA trace. 


and the incremental change in area AA of the oral cavity in response 
to the contained sound pressure. 

Several things are notable. In the sound spectrogram, notice the 
intense initial vowel /a/ with relatively elevated pitch (about 120 Hz) 
and with natural formant transition into the stop. Voicing continues 
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Fig. 7—Behavior of the cord/tract model when the sidewalls are made rigid. The 
input control functions are identical to those of Fig. 6. The synthesized syllable is 
per . Note especially that the cord-oscillator ceases vibration during the mouth 
closure. 


throughout the labial closure, at slightly reduced pitch (about 95 Hz), 
and with the sound output coming solely from the wall radiation. The 
sound level during the lip closure is on the order of 20 dB lower than 
the mouth-radiated vowels. Natural transition into the final vowel 
follows, with voicing at reduced pitch (104 Hz) and intensity. 

The waveforms of the A, and U, oscillations confirm the spectro- 
gram display, as does the waveform of output pressure. The wall- 
radiated sound is dramatized by examining the incremental area 
change in the yielding-wall oral cavity. The area perturbation is seen 
to follow pitch-synchronously the glottal pulses of U,. 
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It is instructive to contrast this soft-wall behavior with that which 
obtains when the tract is made hard-walled; 1.e., by letting Z, >. 
This behavior, for exactly the same input control data, is shown in 
Fig. 7. 

Now, because the tract walls do not yield and permit enlargement, 
the transglottal pressure is rapidly diminished during the labial closure, 
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Fig. 8—Control functions and synthesized sound for the syllable /o’pa/. Sound 
generation for the voiceless consonant is produced from the distributed random noise 
sources that are modulated by the Reynolds number for every network section. 
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and cord oscillation rapidly ceases during the /b/ consonant. Also, no 
sound is radiated from the tract walls, and only silence prevails during 
the lip closure. Offset and onset of cord oscillation, with lip closure and 
release, appears abnormal when compared to transillumination data 
taken on human vocalization. This latter factor may be more important 
perceptually than the actual absence of sound during lip closure. 

Dynamic behavior for a voiceless intervocalic labial stop is displayed 
in Fig. 8. The syllable is /a’pa/, with stressed second vowel. Again, 
control function input is indicated by the top three traces. Only mouth 
area A,, is again displayed, and cord tension Q is held constant. Note 
now, however, the A,. control effects voiced-voiceless switching by 
moving from a value that sustains cord oscillation to one that does not. 

The spectrogram of the sound output shows the low-intensity, low- 
pitch initial vowel with natural formant transitions into the stop. Cord 
oscillation ceases during the closure because the cords are overtly 
pulled apart. (The lateral and posterior crico-arytenoid muscles ac- 
complish this in the human larynx.) The cords come back together 
as the lip closure is released, and oscillation starts with an abrupt 
bounce that is quite characteristically seen in glottal transillumination 
data on humans.!® Natural formant transition is made into the final, 
high-intensity, relatively-higher-pitched vowel. 

The U, flow continues without cord oscillation through the lip 
closure, as the tract wall yields and enlarges the volume forward of 
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Fig. 9—Regions of stable oscillation for the vocal-cord model. 
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the cords. As the lips release, the U, flow reflects a relatively large de 
component before oscillation commences. This flow is the source of 
aspiration in the consonant release. 
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Fig. 10—Behavior of the vocal-cord oscillation as a function of subglottal pressure. 
The dynamics of tract motion and the control of the vocal cord neutral area are the 
same for each diagram. Subglottal pressure for conditions (a) and (b) correspond to 
initial stressed vowel, whereas condition (c) represents final stressed vowel. Note 
especially the delayed voicing onset in (b). 
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The automatic turbulence generation is indicated by the lower trace 
in Fig. 8, which is the squared Reynolds number for the volume flow 
at the lips. As discussed previously, turbulence (noise) intensity is 
monotonely related to this function, in excess of a threshold value.' 
The high spike in R,,, at about ¢ = 0.3 s indicates a turbulent burst 
of noise with approximately this amplitude envelope. The sound out- 
put pressure waveform and the spectrogram show the result of this 
turbulence generation. The result is consistent with aspirated releases 
seen in the /p/ consonant. Furthermore, auditory assessment of the 
synthesized sound indicates a natural-sounding syllable. 

This synthesis also highlights the importance of the A,. control for 
switching between voiced and unvoiced sounds. A more detailed indica- 
tion of this behavior is shown in Fig. 9. Three distinct regions of stable 
cord behavior are indicated. For given cord parameters, stable behavior 
is determined by the interplay of P,; and Ago. 

An additional examination of dynamic behavior dramatizes the so- 
called delayed voicing onset. The syllable /apa/ is generated with the 
A, and A,» controls shown at the bottom of Fig. 10. The cord tension, 
Q, is maintained constant. Lung pressure, P,, however, is varied to 
correspond to initially stressed vowels (conditions a and b) and a 
finally-stressed vowel (condition c). Notice especially in condition b, 
the initially rising, then abruptly falling P, conspires with the first 
opening, and then closing A,, control to produce substantial delay in 
the resumption of voicing. This is found characteristically in human 
speech.!® 


V. AUTOMATIC GENERATION OF CORD/TRACT CONTROL 


Ultimately, we wish to use the cord/tract model as an end-organ 
for speech synthesis. What are the prospects for obtaining the necessary 
controls automatically by rule? 

In recent work on synthesis-by-rule, Coker and Umeda” generated 
synthetic speech from printed text using programmed algorithms for 
articulatory dynamics and for speech prosody. Their speaking machine 
includes a pronouncing dictionary, a syntax and prosody analyzer, 
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Fig. 11—Automatic generation of control functions for the cord/tract synthesizer 
from printed text input. 
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Fig. 12—Example of automatic synthesis of a voiced sentence from printed text using the cord/tract synthesizer. 
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_ Fig. 13—Example of automatic synthesis from printed text for a sentence contain- 
ing voiceless consonants. 


and a dynamic model of vocal-tract shape. The text synthesis program 
calculates several functions that can be transformed into the param- 
eters needed for the control of our cord/tract synthesizer. The sequence 
of conversions is illustrated in Fig. 11. As determined from the Coker- 
Umeda machine, overall sound intensity can be related to P;, voice 
pitch to Q and P,, voiced-unvoiced switching to A,., and tract shape 
to N and A(z). 


TOP VIEW 





VOCAL TRACT 


FRONT VIEW 


LUNGS 


Fig. 14—Format of the computer movie illustrating dynamic behavior of the vocal- 
cord model. 
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(f) 


[eo ey 


Fig. 15—Frame sequence from the computer movie illustrating vocal-cord vibration. 


With the collaboration of Coker and Umeda, we have made an initial 
trial at synthesis of connected speech by making a transformation of 
the prosody and area-function output of their text-synthesis machine. 
An illustration of this first attempt to marry the two systems is shown 
in Figs. 12 and 13. Figure 12 includes plots to show how the machine- 
determined values of voice pitch frequency and intensity are trans- 
formed into the Q and P, parameters required by our synthesizer. The 
spectrogram of Fig. 11 shows completely automatic synthesis of a 
voiced sentence. Figure 12 shows automatic synthesis of a sentence 
containing voiceless sounds. 


VI. SLOW-MOTION COMPUTER PICTURES OF CORD AND TRACT BEHAVIOR 


To aid in visually assessing the complex control and interaction of 
the model components, we programmed high-speed microfilm dis- 
plays of the cord and tract motion. The 16-mm movie film, when shown 
at 24 frames/s, corresponds to a 100:1 slowdown of real time. One 
can, therefore, examine detailed cord motion and _ cord/tract 
interactions. 

One display shows details of the two-mass vocal-cord model under 
dynamic control. The film format is given in Fig. 14 and shows simul- 
taneously a top view of the glottal opening and a front (anterior- 
posterior) view of the two-mass cord model. Some prints of frame 
sequences are given in Fig. 15. The time between displayed frames is 
20 ms. 

A second display, given in Fig. 16, shows a schematized side view 
of the whole vocal system. The vocal tract is simplified to four cyl- 
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LUNGS 


Fig. 16—Format of the computer movie showing dynamic articulatory relations 
between lung pressure, cord motion, and tract shape. 
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Fig. 17—Frame sequence from the computer movie illustrating dynamic behavior 
of the cord/tract synthesizer. The sequence is taken from the synthesis of /a’pa/. 





indrical sections only, but with lengths and areas that change with 
time. The magnitude of subglottal pressure is represented by the 
elliptical contours in the lung volume which expand or contract with 
time. Figure 17 shows a sequence of motion frames, spaced by 20 ms, 
for generation of the syllable /a‘pa/.* 


VIL CONCLUSION 


Initial experiments with this formulation of cord and tract properties 
suggest that the physiologically-based control functions have distinct 
advantages in terms of “built-in” information. That is, much natural 
behavior—such as vagaries of voicing onset and offset, fine-structure 
pitch fluctuations occasioned by tract motion, and voicing behavior dur- 


“The data of Fig. 17 correspond to those of Fig. 8. 
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ing occlusion—is produced automatically in the model. In other words, 
faithful modeling of significant physiological parameters leads to in- 
put control data that can be rather parsimonious. It is therefore not 
necessary to describe input commands with intricate, high-information- 
rate detail. The model is able to generate many of these intricacies of 
natural behavior from relatively simple input control. 

If continued work proves the cord/tract formulation to indeed 
possess the desired physiological constraints and attributes, the 
synthesis approach would also seem promising as a relatively sophisti- 
cated end-organ synthesizer which could be driven by models of 
prosody and articulation, such as provided by the Coker-Umeda text- 
synthesis system. This is an ultimate long-range goal. 

Further than this, however, the model promises some extensive 
potential for studying the dynamics of real speech. Feasibility is 
presently being examined for automatically adapting the model’s 
synthetic output to match real speech waveforms (for example, in a 
least-squares sense). Gradient-climbing adaptive algorithms are being 
examined for this.7 Obvious difficulties are model nonlinearities and 
multiple local-minima traps which may be encountered. Continued 
work will determine whether these analytical questions can be solved. 

Finally, since the present cord/tract synthesis model incorporates 
the technique for automatic generation of turbulence devised earlier,‘ 
this feature permits detailed study of the remarkably delicate articu- 
latory timing the human employs in transitions between voiced and 
voiceless sounds. The cord/tract model therefore fills a critical need 
for a framework within which to organize and assess articulatory 
measurements now being accomplished." 
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The development of the voice-frequency active filters for the D8 channel 
bank is described. These filters are the first single-substrate RC active filters 
using thin-film tantalum RC and silicon integrated-circuit technology to be 
produced on a large scale by Western Electric. To create complete con- 
fidence in both the design and the new technology, an extensive model 
building and testing program was undertaken. In addition, continuous 
interaction with manufacturing engineers resulted in a design that 
facilitated the introduction of this new technology in a large-scale manu- 
facturing environment. 

Significant advances were made in reducing the complexity of tuning 
active filters. In fact, the tuning and testing procedure has been adapted 
for use with a totally automated computer-controlled test set. Furthermore, 
a Monte Carlo statistical simulation of the manufacturing process 
of the filters was developed. This model includes tolerances of the manu- 
factured components, test set errors in measuring gain and component 
values, resistor adjustment accuracies, and temperature and aging be- 
havior of the components. This computer program has been an invaluable 
tool in determining the requirements for tuning, testing, and optimizing 
the final design for minimum manufacturing cost. 


I. INTRODUCTION 


Resistance-capacitance active filters are a relatively new addition to 
the family of frequency selective networks. Although active filters have 
been in existence for 20 years, they have not been widely used because 
passive filters have been less expensive. However, the concurrent de- 
velopment of silicon integrated circuits, tantalum thin-film technology, 
and automated computer-controlled test sets have made Rc active 
filters a practical alternative to passive filters at voice frequencies. 

This paper deals with the design and development of the voice- 
frequency Rc active filters used in the D3 channel bank. These filters 
are realized in thin-film tantalum technology with beam-leaded silicon 
integrated-circuit operational amplifiers on one ceramic substrate.* 
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New analysist:> and modeling techniques, more efficient optimization 
algorithms,’ the concept of statistical design,’ and more flexible 
machine aids to physical design gave the network designers powerful 
tools with which to attack the design process. 


il. DEVELOPMENT PROCESS 


The development process involves a number of interrelated steps. 
Neglecting any of the steps can result in a design that either does not 
meet requirements or is unmanufacturable at a reasonable cost. 
Traditionally, the first step is a careful analysis of the filter require- 
ments. This will determine the order of the filter and, possibly, the 
technology needed for the physical realization (e.g., Lc, Rc active, or 
mechanical technology). Next, the initial design is undertaken. This 
step requires accurate models for the components and general-purpose 
analysis and synthesis routines. Unfortunately, when initial designs 
using new technologies are breadboarded, requirements are very 
rarely met. In addition, the physical realization is either too large 
and/or the network is too costly. Thus, we must optimize the network 
to meet not only electrical requirements but also size, environmental, 
and cost restrictions. Next, sensitivity studies must be undertaken to 
determine both the viability of the design and whether or not manu- 
facturing adjustments (tuning) are necessary. Since the designer 
cannot afford to wait 10 or 20 years to see if his design meets perform- 
ance objectives in the field, a statistical simulation of the manu- 
facturing process and field performance is useful. This simulation can 
analyze the effects of component and adjustment tolerances, simulate 
the tuning procedure and temperature and aging effects, and statis- 
tically optimize the design to give both the greatest manufacturing 
yield and a prediction of the end-of-life performance. 

In the design of filters, early attention to the physical realization 
and tuning and testing for manufacture is important. Lack of atten- 
tion to these details can result in an unmanufacturable design and/or 
a design that has excessive cost. In addition, when new technology is 
introduced, the designer must closely follow the initial manufacture 


of his design to aid the manufacturer in overcoming initial production 
hurdles. 


ill. FILTER ENVIRONMENT AND REQUIREMENTS 


The D1 and D8 channel banks! are the terminals for the Tl pcm 
repeatered line.*.» The T1 pcm line carries 24 telephone conversations 
over two cable pairs. A simplified sketch of the voice frequency end of 
the D3 channel bank is shown in Fig. 1. A voice-frequency analog 
signal from the switching equipment passes through a hybrid trans- 
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Fig. 1—Simplified model of pcm channel bank. 


former. It is then amplified and bandlimited by the transmitting low- 
pass filter. Subsequently, it is converted into a pulse-amplitude modu- 
lated (pAM) signal by a JFET (junction field-effect transistor) switch 
operating at an 8-kHz rate. The pam signals from the 24 channels are 
sequentially encoded into a binary bit stream (1.544 megabits/s) and 
sent out as Pcm signals over a line. On the receiving end, the incoming 
PoM signals are sequentially decoded and converted to PAM signals. 
The analog signal for each individual channel is then recovered from 
the PAM signal by passing it through an interpolating receiving low- 
pass filter. The analog signal then passes through the hybrid trans- 
former to the switching equipment and subsequently to the subscriber. 

The transmitting band-limiting filter must present a good impedance 
to the hybrid transformer and negate the effect of the switched load 
impedance. In addition, it must pass frequencies between 200 and 
3200 Hz with less than +0.09-dB ripple and suppress frequencies 
above 5 kHz by at least 30 dB. At half the sampling rate of 8 kHz 
(or 4 kHz), there must be at least 15 dB of suppression. An Rc active 
filter is ideal for this voice-frequency application, since the operational 
amplifier output will be impervious to the time-varying load, and the 
passband ripple performance is not degraded by inductor losses. The 
above requirements can be met with a fifth-order Cauer-Chebyshev 
filter. To obtain a low-cost, highly reliable filter conducive to high- 
volume manufacture, thin-film tantalum technology*® was chosen for 
the physical realization. 
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One method of realizing a multiple-order active filter is cascading 
noninteracting lower-order sections. This realization consists of two 
stages. The first stage is a modification of a low-pass notch section 
(Fig. 2) developed by J. J. Friend.’ It is a differential-input single 
operational-amplifier section. The second stage is a twin-T notch 
section" that has been modified by the addition of an rc network 
(Fig. 2) to include the pole on the negative real axis. Their respective 
transfer functions are 


Tx(s) = _ tea (1) 


Wyl 

s? + —— § + wey 
Qn 

and 


s? -+- wep 


T2(s) = (+S tenere 


(2) 


The nominal performance of each section and the overall design 
is given in Fig. 3.” 

The smoothing filter on the receiving side of the channel presents 
the designer with a difficult problem. It has a time-varying, rather than 
a time-invariant, generator impedance. This is caused by the JFET 
switch, which operates at an 8-kHz rate (125 us) and is ‘‘on”’ for 3.5 us. 
The switch working in conjunction with the input impedance of the 
network contributes an additional frequency-dependent gain char- 
acteristic to the filtering function of the network. A detailed descrip- 
tion of the design of this filter is given in the next section. 


Ll Hah 
ee Ld 
Fig. 2—Configuration of transmitting filter. 
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Fig. 3—Nominal performance of transmitting filter. 


The design objectives of the receiving filter are as follows :* 


(t) Passband ripple (200 Hz-3200 Hz) < + 0.16 dB. 
(iz) At 3300 Hz: 0.16 dB = gain = — 0.60 dB. 
(17) At 3400 Hz: 0 dB = gain = — 1.20 dB. 
(iv) At 4000 Hz: gain < — 15 dB. 
(v) Above 5000 Hz: gain < — 30 dB. 
(vz) Gain at 1000 Hz: 4.75 dB + 0.02 dB. 


In addition, the filter must absorb the following manufacturing 
and environmental variations: 


(t) Switch “on time’’: 3.1 us to 3.7 us. 
(iz) Switch ‘‘on resistance”: 50 2 to 200 Q. 
(ii) Temperature: 0°C to 60°C. 
(iv) Aging: 20-year life. 
IV. INITIAL DESIGN OF THE RECEIVING FILTER 
The configuration and nominal performance of the receiving filter 
are given in Figs. 4 and 5, respectively.3.> If the switch were not pres- 
ent, the requirements in Section III could be met by a fifth-order 


* All gains and losses are relative to 1000 Hz. 
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Fig. 5—Nominal performance of receiving filter. 


Cauer-Chebyshev transfer function of the form 


T(s) oe 
[s?+ (wp1/Qp1)8 top ILs?+ (wp2/Qp2)S+w5e | (sta) 


The first pole-zero pair can be realized with a second-order tuned 
twin-T.14 The real pole at s = — a could be realized with a first- 
order section, although an operational amplifer can be saved by com- 
bining this pole with the second pole-zero pair to form a third-order 
twin-T section. 

To maintain the gain level and reduce the sensitivity to switch ‘‘on 
time” and “‘on resistance” variations, the input stage must provide 
some “holding” function.’ An ideal sample-and-hold network would 
accomplish this, but it would require an additional operational 
amplifier. Therefore, a holding capacitor was added to the input 
stage. When the switch is on (for a nominal 3.5 us), the capacitor 
C, is charged. During the off time (121.5 us), the capacitor will slowly 
discharge through the network if the input impedance is at a high 
level (>100 kQ). The switch working in conjunction with the high 
input impedance of the first stage acts as a lossy sample-and-hold 
network with a (sin x)/a-type frequency-dependent gain character- 


(3) 
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istic. This gain characteristic has minima at multiples of the 8-kHz 
switch rate (£/E in Fig. 5). 

To design this filter, we could first assume that no switch or holding 
capacitor is present. An initial network would then be developed, to 
which the switch-and-holding capacitor would be added. Because of 
the roll-off introduced by the lossy sample and hold, further optimiza- 
tion would be required. Using a state-variable switched-filter analysis 
program in conjunction with a general-purpose optimization pro- 
gram,*® the desired frequency characteristic is obtained (Tig. 6). 

The program output gives the element values of the first stage of the 
filter (including the holding capacitor) and the transfer function co- 
efficients of the second stage of the filter. To realize the second stage, 
the element values must be described as a function of the transfer 
function coefficients. The transfer function for a tuned, twin-T, second- 
order section with a preceding kc stage can be written as 


3 + w 


Ee) AST Se Giggs enle a) (4) 
or 
_ & + wi 
1) = RT ast + Ap + a . 
If we define 
6 = 1+ (RA/RB), 
a= (2CL/C1), 
cs = (CL/C), (6) 
n= (2R1/RL), 
r= R/RL, 


and assume the twin-T is symmetrical and tuned, the coefficients of 
(5) can be written as 


Ko Bot ara) = 
= oa (8) 
4i= ata | get me (2-6 + 2S) 
+po(?-etatZ)|-S+e © 
4i= Gaye | Wop + om (2-6 +34) 
“+ Rae 27 p+n)|=3+ Sta (10) 
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RC\1l+e+c/ (R1 Cl) 


Since there are more parameters (elements) than constraints (co- 
efficients), the designer can optimize both the absolute values of the 
elements and the element value ratios, in addition to matching the 
required coefficients. An interactive computer program!® was used to 
obtain the element values of the third-order twin-T. 


V. SENSITIVITY CONSIDERATIONS 


An important concern in the design of filters is the degradation of 
the frequency response because of environmental and element value 
changes. Some portions of this degradation can be controlled by speci- 
fying components with tight tolerances, or by carefully tuning the 
filter. Others can be controlled by matching the temperature co- 
efficients of the components used. However, there are some uncontrol- 
lable changes, e.g., aging, adjustment tolerances, and measurement 
tolerances. The sensitivity of the filter design to these changes will 
determine whether a particular technology can be used in the realiza- 
tion, the amount of tuning required, the margin requirements of the 
design, the tolerance of the components, and, subsequently, the cost 
of the filter. 

Many approaches can be taken for analyzing the sensitivity of the 
filter to element changes. One is the degradation of the frequency 
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response caused by all elements taken one at a time. This method is 
very useful in determining which elements require tight tolerances and 
whether tuning or tracking is required. Another is the variation of the 
frequency response under true manufacturing and environmental 
variations. This would involve all components varying according to 
some statistical description. 

As an illustration of the first type of sensitivity, Fig. 7 shows the 
degradation of the frequency response of the receiving filter when each 
capacitor is increased in value by 1 percent. Although the response 
degrades significantly with these capacitor variations, the degrada- 
tions tend to cancel if the capacitors vary in the same direction (track). 
In fact, as can be seen from eqs. (5) through (11), if the second-stage 
capacitors track perfectly, i.e., all capacitors increase by 1 percent, 
the only effect on the frequency response of that stage will be a shift 
in frequency. 

Also, from eqs. (5) through (11), it can be shown that if, in addition, 
the resistors vary in the opposite direction from the capacitors, then 
even this degradation is cancelled. Another point evident from Fig. 7 
is that even +1-percent capacitor variations would cause the response 
to miss the frequency requirements. Since thin-film capacitors can 
only be manufactured to --5-percent tolerances, it is evident that 
initial tuning of the network is necessary. 

Fortunately, thin-film resistors and capacitors have opposite tem- 
perature and aging characteristics. Thus, once tuned, the frequency 
response is quite stable. However, since the tracking is not perfect, 
further analysis must be undertaken.!© This is explained in Section 
VO. 


VI. TUNING 


Many elegant tuning procedures have been developed for hybrid 
integrated circuits.!’— Most rely on gain and/or phase measurements 
at a number of frequencies (one for each transfer function coefficient) 
and adjustment of element values in some algorithmic fashion. Usually, 
some form of descent algorithm, depending upon component sensitivity, 
with a number of iterations is used to solve this multiparameter opti- 
mization problem. Thin-film networks contain critical constraints. The 
capacitors cannot be adjusted and the resistor values can only be in- 
creased either by anodization” or laser trimming.” 

For the above network, nine constraints must be satisfied. They are 
the zero frequencies of each stage, the frequency of the real pole of the 
second stage, the frequencies and q’s of the two complex poles, and the 
1-kHz flat gain. In addition, the manufacturing deviations of the 
switch ‘‘on resistance’ and the holding capacitor must be compen- 
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Fig. 7—Sensitivity of frequency response to capacitor changes of +1 percent. 


sated for. Consultations with manufacturing engineers indicated that 
most iterative approaches that involved measurements of the per- 
formance at a number of frequencies were impractical from a through- 
put viewpoint, i.e., they took too much time. 

From the manufacturing viewpoint, two points are critical: (2) the 
adjustment sequence must be fast and (iz) the adjusted network must 
meet the given frequency requirements. To accomplish both these 
aims, the tuning procedure shown in Fig. 8 was developed. It con- 
sists of two parts. The first is parametric; that is, it depends only upon 
component measurements and component adjustments. The second 
is functional; the gain is measured at one frequency for each stage, and 
an adjustment is made to bring the gain at that frequency to its 
nominal value. Finally, the voltage divider network at the output 
is used to trim the absolute through gain at 1 kHz to a predetermined 
level. 
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Fig. 8—Tuning procedure. 
















ADJUST GAIN 
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To tune the first stage of the filter, the switch-and-holding capacitor 
are initially neglected. The remaining part is a second-order twin-T 
notch filter whose transfer function is given in the appendix. Since the 
capacitors Cl, C2, C3, and CZ cannot be adjusted, their values are 
measured. To obtain a null at the prescribed frequency, eqs. (22), 
(25), and (27) dictate that 


1 _ RI+ R2 nee 
(Ci + OD)R3 ~ RIR2C3”? 
1 
p gees ig 
“2 = PLR? 03 OS (13) 


Since the capacitor values are known, the constraints given by (12) 
and (13) can be satisfied if we pick 


1 
Bi me OS Ca) (14) 
and 
Rl C83 
R3= > op4 G3 (15) 
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To set the pole frequency, eqs. (22) and (82) can be used. Thus, 


2 - Lt (2R1/RL) 
> T+ (L/CS) °* 


Since all the variables but RL of (16) are known, RL can be used to 
set w,. Rearranging (16), 


RL = 


(16) 


2R1 
[1 + (CL/CS) ](wp/w2)? — 1 
At this point, we have calculated the resistor values necessary to set 


the null, null frequency, and pole frequency of (28). To set the pole q, 
eq. (30) can be used. It is repeated here for convenience. 


Wp | ct 2CL 


(17) 


1 
a1 =| R3C1* R2CGS8* Ri c3 Cs 


De 
+ gros |/(+ G8): (18) 


At this point, all the capacitors have been measured and F1, R2, R3, 
and RL have been adjusted. Thus, the only additional adjustment that 
can be made is on 8, which is controlled by negative feedback resistors 
RA and RB. In fact, 


RA 
ae Bp (19) 
Solving (18) for 8, we have 


= 1 1 2CL 1 Wp CL\ |. 
B=1+k3C1 {as | a+ area t az |- 2+ oe) (20) 
To adjust 8 to coincide with that calculated by (20), RA and RB 
are adjusted. 

Because we have neglected the holding capacitor and switch ‘on 
resistance,’ it is necessary to functionally adjust this stage. This 
adjustment is done at 2100 Hz, where the desired nominal gain with 
respect to 1000 Hz is known. Through sensitivity studies, it was found 
that the gain deviation from nominal at 2100 Hz was linearly depen- 
dent upon the ratio of the 8 resistors, RA and RB. Hence, either RA 
or RB is adjusted to complete the tuning of the first stage. 

Since there is interaction between the real pole, the complex pole, 
and the transmission zero, the tuning of the second section is more 
complicated. With C1, C2, and C3 (Fig. 5) measured and the zero 
frequency w, given, the twin-T resistors R1, R2, and R3 can be deter- 
mined. With C and CL measured and a, wy, and gp given, eqs. (9) 
through (11) are a nonlinear set of equations in terms of 8, R, and RL. 
Using the nominal values of these parameters as starting guesses, the 
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new values which set a, w,, and gp to their design values are easily 
found. 

After anodizing the resistors R, Rl, R2, R3, RL, RA, and RB to 
their calculated values, a functional adjustment is performed at 3300 
Hz where the nominal gain value is known. This adjustment mops up 
the previous capacitor measurement and resistor adjustment errors 
and operational-amplifier variations. As in the first stage, the resistor 
ratio of the negative feedback network is touched up in a functional 
adjustment. 

Finally, the gain at 1000 Hz must be adjusted to within +0.02 dB. 
To increase the output level, RG2 is increased. To decrease the level, 
RG1 is increased. 


Vil. TOLERANCE ANALYSIS 


The network designer’s work would be cut significantly if inexpen- 
sive components could be manufactured to their nominal values and 
have no variations with time and temperature. However, manu- 
facturing processes are such that the designer must consider variations 
in the component values and characteristics from sample to sample. 
He could do a worst-case analysis in which he would specify tolerances 
that would ensure that, under all component combinations, the final 
manufactured network would meet specifications. In general, this 
would result in specifying tighter tolerances than needed and neces- 
sarily increase the cost of the network. 

To get a realistic estimate of the performance of the filter as it left 
the manufacturing facility and a good prediction of its field perform- 
ance, a statistical simulation of the manufacturing process and en- 
vironmental behavior was developed. Included in this model were 
manufactured element variations, measurement errors, adjustment 
errors, temperature and aging characteristics of the components, and 
switch timing variations. The variables were the tolerances and dis- 
tributions associated with the above errors, the adjustment procedure, 
and the adjustment frequencies. The figure of merit was the highest 
possible end-of-life yield at the lowest possible cost. Figure 9 is a flow- 
chart for this simulation. 

The simulation proceeds in the following manner. First, a set of 
random numbers is generated for a network. Next, the manufactured 
capacitor and JFET ‘‘on resistance” values are calculated. For the 
thin-film capacitors, the absolute tolerance and the tracking tolerances 
are included. The switch ‘‘on resistance” has a nominal value of 125 
ohms and can vary between 50 and 200 ohms. The next step is a 
simulation of the tuning procedure, where the capacitors are first 
measured, the resistor values calculated, and the resistors adjusted. 
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Fig. 9—Flow chart of tolerance analysis program. 


For this parametric adjustment, the errors in capacitor measurement 
and the anodization errors for the resistors are included. During the 
functional touch-up adjustment, gain measurement errors are simu- 
lated. To evaluate the field performance, we use the temperature and 
aging characteristics of the components. 

Finally the network is analyzed for the final element values that 
are determined by adding the temperature and aging deviations to 
the manufactured element values. Then the network performance is 
compared with the specifications at a number of frequencies, and 
statistics are compiled. If the requirements are not met at any one 
frequency, the network has failed. Next, a new network is picked, and 
we repeat the process. 

The statistical performance of this filter is shown in Figs. 10 and 11. 
Only the minimum and maximum deviations at each frequency for 
1000 sample networks are shown. At each frequency, there is a normal 
distribution of the networks’ performance between the minimum and 
maximum deviations. At room temperature, 25°C, 94.5 percent of the 
samples fell within the design objectives, while at 60°C, 89.3 percent 
fell within. To ensure that the mathematical modeling was correct, 
these results were compared with the measured performance of the 
first 1000 filters manufactured. All those filters fell within the given 
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Fig. 10—Overall statistical performance of filter. 


statistical bands. This type of extensive simulation is justified when we 
have a well-characterized technology, automated adjustment and 
testing, and high-volume manufacturing. 


Vill. PRODUCTION EXPERIENCE 


As previously mentioned, the tolerance analysis program has aided 
in obtaining a practical tuning procedure. If temperature and aging 
variations are ignored, the program predicts the filter yield at the end 
of the manufacturing process. The assigned component tolerances have 
a significant effect on manufacturing costs; thus, a trade-off takes 
place between manufacturing cost and filter yield. Ideally, we would 
like to find the set of component tolerances that minimizes manu- 
facturing cost. This problem has been attacked by Karafin” and Pinel 
and Roberts.?3 However, their work has been restricted to networks 
where: (¢) the performance criteria (i.e., gain or loss constraints at 
given frequencies) is met 100 percent of the time; (iz) tunable elements 
are banned, and (777) there is no correlation between the elements. 

In the case of the D3 filters, it is known that the 100-percent per- 
formance criteria restriction will not produce a minimum cost network. 
In other words, if the yield at final test drops from 100 to 96.7 percent 
by increasing the tolerance on the resistors during parametric adjust- 
ment, the cost drops by more than 3.3 percent. Thus far, cost minimi- 
zation has been attacked in a heuristic manner. Figure 12 shows the 
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Fig. 11—Statistical performance of passband at temperature extreme. 


theoretical manufacturing yields for the transmit filter at 0.1, 0.2, and 
0.5-percent resistor tolerances. Current production experience indi- 
cates that widening the tolerance on certain resistors to 0.3 percent 
results in a lower cost, even though the overall yield decreases. 

One might ask, Why not remove all tolerance restrictions on the 
parametric adjustments and accept any filters as long as they are 
functionally adjusted to meet the frequency requirements at a finite 
set of test points? The answer is quite simple. If all resistors were ad- 
justed at the parametric step to, say, 0.5 percent tolerance, then only 
50 percent of the filters could meet requirements after functional 
adjustment. Since the silicon is bonded after parametric adjustment, 
a much higher final yield is needed to justify the additional investment. 
In addition, as the tolerance is relaxed at the parametric adjustment 
step, more frequency tests are required at final test. 

The tolerance analysis program also helped to provide a temporary 
solution to a capacitor tracking problem. If the ratio of certain capaci- 
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Fig. 12—Statistical performance of transmitting filter with different resistor tolerances. 


tors is not within specified limits, then some resistors may have to 
be anodized excessive amounts—excessive because the resistors will 
then age poorly. However, the tracking requirement on all sets of 
capacitors need not be the same; thus, in a particular case, with the 
tolerance program, it was possible to demonstrate that the tracking 
requirement on a pair of capacitors could be relaxed and therefore 
solve a temporary production problem. 

A tolerance analysis program can be used only if we have good 
estimates for the various tolerances that affect the manufacturing 
process. This was possible for the parameters that influence the fre- 
quency performance, but not possible for the de gain requirement. 
The 1-kHz frequency gain is required to be adjusted to an accuracy 
of +0.02 dB. This stringent requirement was chosen so that no ad- 
justments would be necessary when D3 channel banks were installed 
in the field. The no-adjustment philosophy saves on installation and 
maintenance cost, but it does require critical tuning for the D3 filters.”* 
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IX. PHYSICAL REALIZATION 


In this case, the physical realization uses single-substrate hybrid- 
integrated circuit (HIc) technology.’ The resistors are tantalum nitride 
thin film, and the capacitors have a base electrode of B-tantalum, a 
dielectric of tantalum pentoxide, and a counter electrode of nichrome 
paladium gold. The conductor paths are also nichrome paladium gold, 
while the operational amplifiers and sFET switch are beam-leaded silicon 
integrated circuits. They are all placed on a 33-by-20 mm glazed 
ceramic substrate with a 34-terminal lead frame (Fig. 138). 
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Fig. 183—Hybrid integrated circuit realization of filter. 
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Fig. 14—Electrical layout of nic realization. 


To minimize both the tuning time and the complexity of the test 
fixtures, the substrate layout was constrained. Thus, all measurements 
are made from the sides and all adjustments from the top of the sub- 
strates. To incorporate the tuning procedure of Section VI, the layout 
must provide the capability of measuring all element values. Thus, 
all network nodes must be brought to the edge of the substrate. 

Traditionally, layouts follow the flow of a network. Thus, for a 
network of this complexity, crossovers are generally required, as well 
as break points, to measure the component values. Crossover and 
break closures are normally processed following filter tuning. To 
eliminate breaks and crossovers and separate the measurement and 
adjustment functions, a unique circuit layout was developed. These 
aims were accomplished by judiciously separating common grounds 
and incorporating break points. These redundant points are brought 
to the substrate edge (Fig. 14) and are subsequently connected when 
the substrate is inserted into a printed wiring board. Three terminal 
capacitor measurements and capacitor electrode sharing were neces- 
sary to make these measurements. To realize this layout with its 
subsequent 11 mask levels required considerable reliance on computer- 
aided graphics.?*.?6 
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APPENDIX 
Transfer Function of Second-Order Twin-T Section 
A.1 Untuned case 
The transfer function of the untuned, unsymmetrical second-order 
twin-T is 
(s + 72)8? + (s + 71)0% 
EGF WD) FOE) a 


with the following definitions: 
CP = Cl + C2, CS = (C1 C2)/CP, 


T(s) = 


RS = R1 + R2, RP = (R1 R2)/RS, 
71 = 1/(R3 CP), tT. = 1/(RP C3), (22) 
5, = T2 — 71, B=1-+ (RA/RB), 
c1 = CL/CS, ri = RS/RL. 
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The polynomials D(s) and E(s) of (21) are 


1— 8 
R3 C1 


s 


RL aE GE | (1 + 1) 


1l+ri 2 
+ (FE) us 


D(s) = e+| peat 4 civaeee 


R2 BOB 





B®) =|ris8+ pets |> 


l+e R3 Ci(1 + C1) 
and 
> 1 
Os = Fa no va wo 
7 *R1 R2C3 CS’ 
_ 6 
k= 1 + C1 
A.2 Tuned case 
When the twin-T is tuned, 
1 = 72, 5, = 0. 
Thus, (21) reduces to 
_ s+ we 
EG) Hs + ais + ae 
Ss? + we 


ie 8+ (wWp/qp)8 + we 


If, in addition, a symmetric twin-T is picked, 1.e., 











Rl = R2 = 2R3 
Cl = C2 = €3/2 
CL = (C1 a) /2 
L= (2R1)/r1, 
then 
T1 + C1 1 
pie pai (2-8 + 24 eure 
a2 = G a rt) ot 
2 1 + C1 2) 
and 





= 1l+n a 
op = (Fy )o 

5 Orr Oey. 
Ye ~ 72 —6+ (n+ c)/2] 
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FIR Digital Filter Banks for Speech Analysis 
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In using filter banks for processing speech signals, it is often important 
that the sum of the individual frequency responses of the bandpass filters 
(composite response) be flat with linear phase. This paper presents a 
technique for achieving flat composite response using linear-phase FIR 
digital filters. The design method is based on some special properties of 
FIR filters designed by the windowing method. Excellent response char- 
acteristics can be achieved with complete flexibility in choosing the center 
frequencies and bandwidths of the individual filters. 


I. INTRODUCTION 


Filter banks are used to perform short-time spectrum analysis in a 
variety of speech processing systems.'~* Typically, a set of bandpass 
filters is designed so that a desired portion of the speech band is 
entirely covered by the combined passbands of the filters composing 
the filter bank. The outputs of the bandpass filters therefore are con- 
sidered to be a time-varying spectrum representation of the speech 
signal. If special care is taken in the design of the bandpass filters, it is 
possible to reconstruct a very good approximation to the input speech 
by simply adding together the outputs of the bandpass filters.® This 
is the basic principle of a variety of vocoder systems. 

Since the bandpass filters are linear systems, we can characterize the 
behavior of such filter banks by considering the composite frequency 
response when all the outputs are added together. Since, ideally, the 
output should be equal to the input, then we desire that the composite 
frequency response have constant magnitude and linear phase in the 
desired band of frequencies. This criterion, together with specifications 
on the desired bandwidths of the individual frequency channels, forms 
a meaningful basis for the design of filter banks for speech analysis. 

An earlier paper® showed that careful attention to the relative phases 
between channels is important in achieving a flat composite frequency 
response. That paper, which was concerned primarily with filter banks 
composed of infinite impulse response (11R) digital filters, described a 
method of obtaining flat composite frequency response by a relatively 
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simple adjustment of the relative phases of the channels. This method 
was later applied to the design of a speech analysis/synthesis scheme 
in which finite impulse response (FIR) digital filters were used.’ Using 
this method, excellent overall response can be obtained for both Rr 
and FIR digital filters in filter banks in which the center frequencies are 
uniformly spaced. However, the method is not easily extended to 
nonuniformly spaced filter banks. 

In the present paper, we describe a different approach that is not 
limited to the design of uniformly spaced filter banks. The method 
exploits some special properties of linear-phase rir filters and thus 
cannot be applied very successfully to the design of 11R filter banks. 
We first discuss the basic design principles, and then show some design 
examples. We conclude with a discussion of some computational 
considerations of Fir digital filter banks. 


il. DESIGN METHOD 


FIR digital filters are attractive for design of speech filter banks for 
several reasons. First, such filters can be designed to have precisely 
linear phase simply by imposing the constraint 


h(n) = h(N — 1 — n) OsnSN-1 (1) 


(on each individual filter band*), where h(n) is the impulse re- 
sponse of the filter and N is its length in samples. This means that the 
criterion of linear phase for the composite filter bank response is 
trivially met if the individual filters have identical linear-phase 
characteristics. Therefore, it is possible to focus attention on achieving 
arbitrary frequency selective properties for the individual filters and on 
obtaining the desired flat response for the composite filter bank. The 
second great advantage of Fir filters is that a variety of design methods 
exist ranging from the straightforward windowing method®.’ to itera- 
tive approximation methods that allow great flexibility in realizing 
complicated design specifications.® 


2.1 FIR bandpass filters 


The bandpass filters that we shall consider have impulse responses 
of the form 


hi(n) = hiz(n) cos (wexnT’) OsnSN-1 
= 0, otherwise, (2) 


where hj;,(n) is the impulse response of the kth linear-phase low-pass 


_ “It is assumed, for simplicity, that the impulse response of each bandpass filter 
is of duration N samples, although it is trivial to remove this restriction by adding 
appropriate delays for each channel. 
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Fig. 1—Implementation of a typical bandpass channel. 


filter. This particular form for the impulse response is motivated by the 
fact that, in some vocoder applications,?* each bandpass channel is 
implemented as shown in Fig. 1. The overall impulse response of the 
system of Fig. 1 from input x(n) to output y(n) is easily shown to be 
given by eq. (2). 

The spacing of the individual channels of the filter bank is deter- 
mined by the choice of the set of center frequencies, w,,, which is in 
turn determined by the desired frequency resolution of the filter 
bank. The frequency selectivity of each channel is determined by the 
frequency response characteristics of the prototype low-pass filters 
hu,.(n). Since phase considerations can be simply avoided by designing 
all the bandpass filters to have the same linear phase, we can focus 
our attention entirely on designing a set of prototype low-pass filters 
that have the desired individual frequency selective properties and 
that give the flattest amplitude response for the composite set of 
bandpass filters. 


2.2 Low-pass filter design 
The window design method appears to have a number of advantages 
for design of the prototype low-pass Fir filters. This method is de- 
picted in Fig. 2. First, a desired ideal low-pass filter of the form 
Hu;,(e7) = e—ien? | cw | < W pk 


0, otherwise, (3) 


is defined by choosing the cutoff frequency w,,. Note that, for sim- 
plicity, we have omitted in the figure the linear phase term 
exp (— jwnol’) corresponding to a delay of n> samples. The value of 
Mo required is no = (N — 1)/2. This means that, if N is even, the 
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Fig. 2—Windowing technique for a low-pass design. 


delay corresponds to a noninteger number of samples. The ideal im- 
pulse response for the kth channel is, therefore, 


_ 1 f* ar jjont a, iD Looe(nT — mo&P)) 
hax(n) = on i eno? ejonTdyy = ane es mae =m) (4) 


Of course, this impulse response is infinite in extent and must be 
truncated to obtain an Fir filter. This is done by defining 

hu(n) = w(n — no)har(n), (5) 
where w(n) is a window function and hi,(n) is the impulse response of 
the kth prototype low-pass filter. The length of the window, denoted 
by N, can be either an even integer (VN = 2M) or an odd integer 
(N = 2M + 1). Figure 2 shows the case when N is odd. 

The result of multiplying the ideal low-pass impulse response by the 
window corresponds to a convolution in the frequency domain of the 
ideal frequency response and the Fourier transform, W(e%*"), of the 
window; i.e., 


x{[T 
H,(e#7) = = i; Has,(e77) W (e4e-©7) da. (6) 
20 J_x/T 
The result of this convolution is depicted in Fig. 2. It can be seen that 


534 THE BELL SYSTEM TECHNICAL JOURNAL, MARCH 1975 


the main effects are the introduction of a smooth transition between 
the passband and the stopband and the introduction of ripples in the 
passband and stopband regions. The properties of this approximation 
are depicted in Fig. 3. If w, is larger than the width of the ‘‘main lobe”’ 
of W(e”7), then the following set of properties are generally true: 


(t) The transition region, Aw, is inversely proportional to N. 
(it) The function H(e’*T) is very nearly antisymmetric about the 
point (wy, 0.5). 
(iit) The peak approximation errors in the passband and stopband 
are very nearly equal. 
(tv) The approximation error is greatest in the vicinity of w,, and 
it decreases for values of w away from w,». 


The above properties of the windowing design method are true of 
all the commonly used windows. However, Kaiser has proposed a 
family of window functions that are very flexible and nearly optimum 
for filter design purposes.® Specifically, the Kaiser window is 


Iola vl — (n/no)? | 
I o(a) 

= 0, otherwise, (7) 

where no = (N — 1)/2 and J,[.- ] is the modified zeroth-order Bessel 

function of the first kind. By adjusting the parameter a, one can trade 


off between transition width and peak approximation error. Further- 
more, Kaiser’ has formalized the window design procedure by giving 


w(n) = || <= No 


H (ei @T) 


0.5 


: a (A — rf 


G 
md 


Fig. 3—Resulting low-pass design from windowing. 
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the empirical design formula 


az —20 logio 6 — 7.95 ats i. 


Mes 14.36Af (8a) 


where N is the filter order, 6 is the peak approximation error, and Af 
is the normalized transition width 


ASS 5a (8b) 


To use this formula, we fix 6 and Af at values that provide the desired 
frequency selectivity. Then eq. (8a) can be used to compute N, and 
the parameter a can be computed from the equation’ 


a = 0.1102(—20 logis 6 — 8.7), — 20 logio 6 > 50 
= 0.5842 (— 20 logio 6— 21)°-4 
+ 0.07886 (— 20 logio 6 — 21), 21 < — 20 login 6 < 50. (9) 


In the present application of this design method, the choice of 6 and 
Af depends upon the specifications of the bandpass filters that con- 
stitute the filter bank. 


2.3 Filter bank design 


To design a filter bank using Fir filters, we must first determine the 
range of frequencies to be covered by the composite response. Let us 
assume that these are denoted wmin and wmax, Where wmax S a/T. 
Now, if there are a total of N, filters, we must choose the bandwidths 
and center frequencies so that the entire range of frequencies 
®min S w S wWmax 18 covered. This is depicted in Fig. 4 for the case 


N; = 8. This figure shows the ideal responses for each bandpass filter ; 


[Hak (eJ@T)| 





OmMin QMAX 


Fig. 4—A typical nonuniform filter bank. 
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i.e., as would be obtained if windowing were not required. In general, 
it is clear that 


Ny 
@max — Wmin = 2, 20ps (10) 
=1 
and 
k—1 
Wek = Wmin + ys 2wWpm + Opk k = 2 
m=1 


= eben k= 1. (11) 


If all the filters have the same bandwidth, i.e., w,, = wo, then it is 
easily seen that 
@max ~~ Wmin— (12) 


2N; 


Wo = 


Alternatively, if the bandwidths are to increase exponentially; e.g., 
Wk = 2* Iwo, then 


Wmax — Wmin (13) 


2(257 = 1) 


Wo = 


The center frequencies can be found in either case by using eq. (11). 
The choice of peak approximation error depends upon how much 
stopband attenuation is deemed necessary in a given application. 
Typical values of —20 logis 6 would most likely be between 40 and 
60 dB. Using eq. (9), the appropriate value of a can be computed. 
Finally, the normalized transition width Aj must be fixed to compute 
N from eq. (8a). Again, the choice of Aw (or Af) is governed by con- 
sideration of the desired frequency selectivity for the individual filters. 
Clearly, the transition width Aw; should not be more than 2wx. 

In the filter bank context, we shall require that Aw be the same for 
all filters so that we can take advantage of property (77) of Section 2.2. 
That is, if all the filters have identical transition regions and, further- 
more, if these transitions are antisymmetric about the crossover 
points, then we can expect that the sum of the frequency responses 
will be very close to unity. This is illustrated in Section III. 


lil. DESIGN EXAMPLES 

In this section, we illustrate the use of the principles established in 
Section II with examples of both uniform and nonuniform filter 
banks. For all the examples, the sampling rate is assumed to be 9.6 kHz. 
Example 1 


Suppose that we wish to design a bank of 15 equally spaced filters 
that covers the range 200 to 3200 Hz. Then, using eq. (12), we find 
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that the cutoff frequency for all the low-pass filters* is 


fo= _ = 100 Hz. 


Using eq. (11), the center frequencies are 


Wek 


for = 5% = 100k +1) Hz k= 1 





ee 


If we assume that 60-dB attenuation is required outside the transi- 
tion regions of each channel, we find from eq. (9) that a = 5.65326. 
Since the cutoff frequency is 100 Hz for all the prototype low-pass 
filters, the widest transition band that is reasonable is 200 Hz. Using 
this value and —20 logio 6 = 60 in eq. (8a), we obtain N = 175 as 
the lowest reasonable value for N. Note that, if lower attenuation is 
acceptable, then N can be smaller for the same Af. 

The filter bank designed with the above parameters is shown in 
Fig. 5. Figure 5a shows the individual bandpass filters. Note how the 
fall-off in the upper transition band of a given filter complements the 
ascent of the next filter. Also note that adjacent channels cross at an 
amplitude value of 0.5. Figure 5b shows the composite response of the 
filter bank. It is clear that the filters merge together very well at the 
edges of the frequency bands. Indeed, the deviation from unity is less 
than or equal to the peak approximation error, 6 = 0.001, that was 
used in designing the prototype low-pass filters. 


Example 2 


A nonuniform spacing of the filters is often used to exploit the ear’s 
decreasing frequency resolution with increasing frequency. Suppose 
that we wish to cover the same range 200 to 3200 Hz as in Example 1, 
but we wish to use only four octave band filters. That is, each succes- 
sive filter will have a bandwidth twice the bandwidth of the previous 
filter. Using eq. (13), we find that the lowest frequency channel has 
cutoff frequency 

wo 3200 — 200 


Ler pe ore 


In general, the cutoff frequencies of the prototype low-pass filters are 


fon = Gr = fy k= iF 2, 3, 4, 


* For the actual low-pass filter, the response will be approximately 0.5 at w = wp, 
the cutoff frequency of the ideal low-pass filter. 
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INDIVIDUAL FILTER RESPONSES (N = 175) 
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Fig. 5—Individual and composite frequency responses of a bank of 15 uniform 
bandpass filters for N = 175. 


or the bandwidths of the bandpass filters are 200, 400, 800, and 1600 
Hz, respectively. The center frequencies are found, from eq. (11), to 
be 300, 600, 1200, and 2400 Hz, respectively. Again requiring 60-dB 
attenuation, we note that the narrowest bandwidth is 100 Hz, so that 
the smallest reasonable transition width is 200 Hz. This leads again to 
a@ minimum value of N = 175. The filter bank corresponding to these 
design parameters is shown in Fig. 6. In Fig. 6a, again note the rela- 
tionship between the ascending and descending transitions between 
adjacent filters. Particularly note that, since N and a are the same for 
each of the prototype low-pass designs, the shape of the curves in the 
transition region is independent of the bandwidth. Figure 6b shows 
the composite response where the deviation from unity is again less 
than 0.001. 
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INDIVIDUAL FILTER RESPONSES (N =175) 





uJ 
Q 
= 
E 
- 
= 
< 0.5 
0 
0 1000 2000 3000 
FREQUENCY IN Hz 
COMPOSITE RESPONSE 
u; 1.000 
ja) 
— 
a 
a 
a. 
= 
<< 
0.995 





0 1000 2000 3000 
FREQUENCY IN Hz 


Fig. 6—Individual and composite frequency responses of a bank of 4 nonuniform 
bandpass filters for N = 175. 


It is interesting to note that the composite frequency response of the 
filter bank is independent of the number and distribution of the in- 
dividual filters, so long as the same window is used to design all the 
individual filters in the bank. This result can be verified by writing the 
overall frequency response of the filter bank, H(e’7), as 


Nt 
H(e"?) = 2, Hi(e*"), (14) 


which, from eq. (6), can be written as 


H(t) = yi Z 


ko1 27 


w/T 
fo Hane?) W (eV a. (15) 
—7r/T 
Interchanging the order of summation and integration, eq. (15) 


540 THE BELL SYSTEM TECHNICAL JOURNAL, MARCH 1975 


can be written as: 


: T x{[T Ny : . 
He?) = = > Hale?) | W (eio-9)T) dg (16) 
WT Jox/T L k=1 
= W(e#T)®H 7(e%7), (17) 
where 

Ny 
Hy(e??) = > Ha.(e%7). (18) 

k=1 


Equations (17) and (18) show that the overall frequency response of 
the filter bank is the circular convolution of the frequency response of 
the window with the frequency response of the combined ideal band- 
pass filters. As seen in Fig. 4, the combined ideal frequency response 
of the bandpass filters is an ideal bandpass filter from w = wmin to 
® = Wmax, Independent of the number and distribution of the in- 
dividual filters. Thus, the composite filter bank frequency responses 
for the examples in Figs. 5 and 6 are identical because the same 
window was used in both cases and the filters spanned the identical 
frequency ranges. 


Example 3 


Suppose that all the parameters remain the same as in Example 2 
except that we require narrower transition regions. This means that 
a larger value of N is required. In fact, Eq. (8a) shows that N and Af 
are roughly inversely proportional. Figure 7 shows the filter bands cor- 
responding to the parameters of Example 2 except that N = 301 and 
Af = 0.012082 (transition width is 116 Hz). The sharper transitions 
are apparent in Fig. 7a, and Fig. 7b shows that the composite response 
remains very flat. 


Example 4 


We have assumed throughout that the transition width was less than 
twice the smallest low-pass cutoff frequency. In our examples, this 
constraint required that N be at least 175. The result of reducing NV 
below this value is illustrated in Fig. 8. In this case, all the parameters 
were the same as in Examples 2 and 38, except in the case of N = 101 
and Af = 0.0362465. The transition width is 348 Hz, which is much 
greater than twice the cutoff frequency of the first low-pass filter. 
This is clearly in evidence in Fig. 8a. It is clear that reasonable filters 
are obtained for the wider bandwidth filters; however, the lowest 
filter does not attain unity response anywhere in its passband. 

The preceding examples make it abundantly clear that, for suffici- 
ently long impulse responses, the composite filter-bank response can 
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INDIVIDUAL FILTER RESPONSES (N = 301) 
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Fig. 7—Individual and As frequency responses of a bank of 4 nonuniform 
bandpass filters for N = 


be very flat. In Ref. 5, where design techniques for rR filter banks were 
discussed, the best results achieved for the composite response were 
approximately 1-dB peak-to-peak ripple for uniform bandwidths and 
about 2.5-dB peak-to-peak ripple for nonuniform bandwidths. This is 
in contrast to the results of the examples of this section, where the 
peak-to-peak ripple in the composite response was about 0.0274 dB 
for all the filter banks independent of how the bandwidths were chosen. 
This, together with the precise linear phase that is easily achieved, 
makes the rir filter banks superior to what can be achieved for 1R 
filter banks. The price that is paid for this is that rather large values of 
N are required to achieve sharp transitions. However, the values of NV 
used in the previous examples are certainly not unreasonable if the 
filters are implemented by FrtT convolution methods or in special- 
purpose hardware. 
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INDIVIDUAL FILTER RESPONSES (N = 101) 
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Fig. 8—Individual and composite frequency responses of a bank of 4 nonuniform 
bandpass filters for N = 101. 


IV. SUMMARY 


We have discussed a design method for filter banks composed of 
FIR digital filters. The method exploits the linear-phase properties 
obtainable for such filters, as well as the symmetry of the transition 
region that results from the windowing method of design. We sum- 
marized this method of design for the Kaiser window and illustrated 
the filter-bank design method with several examples. These examples 
show that the proposed design method has a great deal of flexibility 
and that excellent response characteristics can be achieved. 
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Design and Evaluation of Shifted-Companion- 
Form Active Filters 
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Two techniques for designing a class of low-sensitivity, follow-the- 
leader, feedback-type active filters have been introduced by Hurtig and 
Laker-Ghausi. The FLF configuration consists of a cascade of second- 
and/or first-order sections, with feedback from each section back to the 
first. This paper presents an approach for designing FLF-type realization 
for all classes of filter functions. The technique is based on a shifted-com- 
panion form of the associated-state equations. Some salient features of 
Hurtig’s primary resonator block, Laker-Ghausv’s follow-the-leader feed- 
back, and the shifted-companion-form techniques are presented below. 


(7) Hurtig’s PRB realizes any all-pole (no finite transmission zeros) 
jilter function. This includes the low-pass, high-pass, and sym- 
metrical bandpass filters without finite zeros. Explicit design 
equations are available, and the individual sections in the array 
are identical. 

(12) Laker-Ghaust’s FLF realizes any symmetrical (including finite 
transmission zeros) bandpass filter function. The sections are 
not constrained to be identical, which allows optimization using 
this degree of freedom. Finite zeros are realized by a summation 
technique. 

(11t) The SCF realizes all types of filter functions, 1.e., low-pass, 
high-pass, bandpass, all-pass, or band-reject filters. Explicit 
design equations are available. The first section can differ from the 
rest, thus allowing some optimization with standardization. Feed- 
forward as well as summation techniques can be used to realize 
the finite zeros. 


Two bandpass design examples using SCF, PRB, and/or Laker-Ghaust 
FLF techniques are given and compared with the low-sensitinty coupled 
(leapfrog) biquad, the conventional cascade biquad, and the passive 
ladder filter designs. The comparison shows that the passive filter gives 
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the best performance with respect to sensitivity to element deviations. All 
the coupled designs are significantly better than the cascade design in the 
passband, with the coupled biquad (leapfrog) design the most signif- 
icantly better. In the stopband, cascade and coupled designs perform 
roughly the same. 


1. INTRODUCTION 


Recently, Hurtig!? introduced a low-sensitivity, multiple-loop-feed- 
back active rc filter configuration for the realization of greater-than- 
second-order voltage transfer functions. The configuration has been 
found to exhibit greatly improved stability over cascaded designs. 
For symmetrical bandpass filters, Hurtig’s structure [called the pri- 
mary resonator block (PRB) configuration] consists of a cascade of 
identical biquadratic bandpass sections (i.e., same pole-frequency and 
pole-Q) with feedback from each section (except the first) back to the 
first section. More recently, Laker and Ghausi?* extended Hurtig’s 
configuration to include symmetrical bandpass filters with finite 
transmission zeros, e.g., elliptic-type filters. In Laker and Ghausi’s 
approach [called the follow-the-leader feedback (FLF) technique ], 
different pole-Q values can be allowed for the biquadratic bandpass 
sections. 

For the prs technique, Hurtig has given a set of explicit equations 
expressing the biquadratic bandpass transfer function and the feed- 
back factors in terms of the coefficients of the all-pole prototype low- 
pass transfer function.? In the FLF approach, Laker and Ghausi used 
a coefficient-matching technique. Because of the nonuniqueness of 
solutions in the FLF approach, Laker and Ghausi further proposed a 
method of choosing the pole-Q values for an optimum design. 

In this paper, we present yet another approach based on a shifted- 
companion form of state variable representation of the voltage transfer 
function for the design of symmetrical bandpass and band-reject 
filters with this structure. In the bandpass case, using the proposed 
method, each biquadratic bandpass section in the cascaded array must 
be identical, with the possible exception of the first. Hence, it in- 
cludes the Hurtig prs configuration as a special case, but does not 
encompass the Laker-Ghausi cases having three or more different 
values of pole-Q. As in Laker-Ghausi’s approach, the design of sym- 
metrical bandpass filters with finite transmission zeros is included in 
the discussion of the shifted-companion form. Similarly to Hurtig’s 
approach, the shifted-companion form also gives explicit design 
formulas as opposed to the coefficient-matching technique used by 
Laker and Ghausi. Furthermore, in the shifted-companion-form 
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design, different realizations (with the same configuration) can be 
obtained by varying the value of a shift parameter. The standard 
companion-form representation® corresponds to the case in which the 
value of this shift parameter is equal to zero. 

In the next section, the shifted-companion-form representation of a 
voltage transfer function is presented. A brief discussion on the optimal 
choice of the shift parameter based on our design experience is given 
in Section III. Two design examples, a three-section Butterworth 
bandpass filter and a three-section elliptic bandpass filter, are given in 
Section IV. The section also compares the sensitivity performance, in 
a Monte-Carlo sense, of the shifted-companion-form designs to the 
cascade biquad and the coupled biquad®” as well as to the passive 
designs. 


Il. SHIFTED-COMPANION-FORM REPRESENTATION OF VOLTAGE 

TRANSFER FUNCTION 

The design technique for the proposed shifted-companion-form 
representation of a voltage transfer function is obtained as follows. 
First, a shift is introduced to the complex frequency variable by adding 
a variable constant a (shift parameter) to the complex frequency 
variable. Second, the resulting shifted-transfer function is represented 
by the standard companion form’ and its corresponding block diagram, 
which has the desired structure. Third, an inverse shift operation is 
made on the standard companion form to determine the proper values 
for the parameters of the structure. 


2.1 Representation of voltage transfer function by a shifted-companion form 
Let the voltage transfer function be given by 


Vout jel 1 Nm-1p™ | t+ Mp + No 


Vin pt d,p +. +dpth |” 


for m<n. (1) 


Vag 


Let the following shifting be made in the complex frequency variable 
of (1): 

p =s-— a, (2) 
where a, the shift parameter, is a real number. Substitution of (2) into 
(1) results in the following shifted-transfer function (see Appendix A): 


3 


Dm—is™—* 
Vos a ve 
yo (s) = "+4, (3) 
i gs” + > An—js"—4 
fai 
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ee (n — k)! 7 
f= BOM Ga iia= pi 


(3a) 


i (m — k)! : 
<— _~ Oh NT EEE Ss =) eek, 
= 2 Ge iano 
7=0,1,-:+-,m 
Note that d, = 1. Alternatively, the a’s and b’s can also be obtained 
by the following implicitly recursive formula: 


a | 
j (n — k)! ; 
do = See a an: eae = 1,2, -°:, 
i= 2, G—hin— pt J _ (3b) 
t = ! 
Nm—i = mee £=0,1, +++, 


irk 
vn G—bim-pit or 


It is well known that a voltage transfer function (with degree n) can 
be represented by a set of state equations in the (standard) com- 
panion form,? Le., 
x = Ax + Doin 
Vout = CX + dvin, 


(4) 


where the state matrix A is of dimension n X n. In the case of eq. (3), 
we have 
x= (1, 2 Tn)' 


—~Gn-1 ~An-2 —GAn-3 *** —1 —~Ao Bi 
1 0 0 -- 0 0 b= | & 
A= 0 1 0 -. 0 Of, (5) 
: Bn 
0 0 0 a | 0 


c = [1, v2, °°°, Yn] 


There are two special cases for eq. (5), A and B. 


Case A: Transmission-zero forming by an input feed-forward tech- 
nique: 

c=[00--- 1] 

b = [6:1 B2 «++ Bn]! 
and 


i—1 
Bnzi-i = Oni — yy Qn-itnyi-j, t= 1,2, °°, 7. (5a) 
j= 


548 THE BELL SYSTEM TECHNICAL JOURNAL, MARCH 1975 


Case B: Transmission-zero forming by summation-of-state-variables 
technique: 
b=[10.--- 0] 


c= [bn—-1 Dn—2 eee bo |. (5b) 


To obtain the shifted-companion-form representation of the voltage 
transfer function of eq. (1), the inverse shift operation, i.e.,s = p +a 
is applied to eq. (4). This is equivalent to the following operation: 


time domain frequency domain ime domain 


inverse shift 


s=pta 


where I is the n X n identity matrix. Hence, a shifted-companion-form 
representation of eq. (1) is* 


7 = A’ + Duin 
eae (7) 
Vout = CY of dvin, 
where 
—GQn-1 — @ ~Gn-2 —Gn-3 *°** 41 —Aq 
1 —a 0 --» Q 0 
A’=A-—al = 0 1 —a --» 0 , (7a) 
0 0 O° we Tk ee 


and the vectors b and c are as given in eqs. (5a) or (5b). 

At this point, it is desirable to change the relative level of the state 
vector y to obtain more convenient values for the gain (i.e., close to 
unity) of the individual biquadratic sections. Mathematically, we let 


y = Kx, (8) 


where K is a nonsingular diagonal matrix. It has been found convenient 
to choose K to have the following form: 


K = diag [a* a? --- @ 1). (9) 


Substituting (8) and (9) into (7) and (7a), the following shifted- 
companion-form representation of eq. (1) is obtained: 


x = Ax + bein 


10 
Vout = cx dvin, ( ) 


“The state vector is changed from x to y. 
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where 














—a a= = An—2 = An-3 fe a1 ao 
n—1 Q A nt an a? 
ee ae a —a 0 ee 0 0 
peri ein 0 a —a oe 0 0 
0 0 0 a —a 
(11) 
b = Kb = | ses Be J Bat 2 | 
Qa a Qa 
é=cK=[0 0 --- 0 1] (11a) 
or 
‘ 1 t 
b= Kb =| O --- O 0| 
a 
¢ = cK = La" bn_1 ab,» oss ab; bo |. (11b) 


Equations (11a) and (11b) correspond to the cases where the trans- 
mission zeros are formed by the input feed-forward and the summation 
techniques, respectively. 


2.2 Block diagram representation of the shifted-companion form 


Transforming eqs. (10), (11), and (11a) into the frequency domain, 
the following set of transfer functions representing the shifted-com- 
panion form is obtained. 





2 u Sh ei y (yy 4 By, 
K0) = eee |- BHO +r | 
Xi(p) = ata | eX) + oe Vial) | > (12) 
for 4 = 2,3, ---,n 


Vout (p) _ X n(p) ae dV in(p) 
Similarly, by transforming eqs. (10), (11), and (11b), we have 


n 


1 Qn-j 1 
Xi(p) = eG ee) | - » = X;(p) + san into) | 


j=2 





1 
pte 


Vout(p) = = ab, _:Xi(p) + dV in(p) 


Xi(p) = 





[eX:s(p)] for i= 2,3,---,n - (18) 


Equations (12) and (13) are shown in block diagram form in Figs. la 
and 1b, respectively. 
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Fig. 1—Shifted-companion form. (a) Feed-forward zero-forming technique. (b) Summation zero-forming technique. 


2.3 Block diagram representation of symmetrical bandpass filters via the 
shifted-companion form 

For geometrically symmetrical bandpass filters, eq. (1), eq. (12) or 
(13), and Fig. 1a or 1b can be taken as the transfer function, the shifted- 
companion-form representation, and the block-diagram representation 
of the corresponding low-pass prototype, respectively. To obtain the 
block diagram representation of the symmetrical bandpass filter 
transfer function, we can apply the well-known low-pass to bandpass 

transformation: ‘ 
_ & + a 


Bs oa 


where 


p= complex frequency for the normalized low-pass function 


s = complex frequency for the actual bandpass function 
wo = center frequency of the bandpass filter (in radians/s) 
B = passband bandwidth of the bandpass filter (in radians/s) 


to Figs. la and 1b. The resulting block diagram representations are 
shown in Figs. 2a and 2b. 


2.4 Block diagram representation of symmetrical band-reject filters via the 
shifted-companion form 

To obtain the block diagram representation of the symmetrical 
band-reject filter transfer function, similarly to the development of 
Section 2.3, we can first apply a low-pass to high-pass transformation, 
then follow with the usual low-pass to bandpass transformation, eq. 
(14). Specifically, this results in the following transformation to Figs. 
la and ib: 


1 & + wf 


1 
pa a s+ (B/a)s +a’ do) 
where 
p = complex frequency for the normalized low-pass function 
s = complex frequency for the actual band-reject function 
wo = Vwew1 (in radians/s) 
B = w. — o (in radians/s) 


w/w: = the lower/upper passband edge frequencies of the band-reject 
filter. 


The resulting block diagram representation can also be shown as in 
Fig. 2, except that 


2 2 
Ti(s) = : oie 


Qn-1 ta s+ [Bs/(ar-1 + a) ] + oF 
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Bs Bs : (b) 


WHERE Ty,(s) = Ti(s) i=2,3,...,n 


s2+B(@n—1ta)st w o2 s2+Barst@ 92 


Fig. 2—Symmetrical bandpass filters via the shifted-companion form. (a) Feed-forward zero-forming technique. (b) Summation zero- 
forming technique. 


and 
St 0 


"2 + (BJa)s + we’ 


Ti(s) = + t= 2,3, +++, 1. 

In passing, we note that Fig. 2 can also represent the symmetrical 
band-reject filter provided the parameters in Fig. 2 are determined 
by making eq. (1) the transfer function of the band-reject’s correspond- 
ing high-pass prototype. 


Il. OPTIMAL CHOICE OF THE SHIFT PARAMETER, a 


For symmetrical bandpass filters (Fig. 2), it is seen that all the 
biquadratic sections, with the possible exception of the first section, 
have a pole-Q value equal to wo/Ba. The pole-Q value for the first 
section is 

Wo 0 


Ba@.¢0) © ‘Bian =-G@— Del 


The value of these Q’s versus a is illustrated in Fig. 3. 

Before we proceed with a discussion on the optimal choice of a, 
two special cases are pointed out. The first is the standard companion 
form which corresponds to the case where a = 0. From Eq. (8a), 


Qn-1 = An-1 —- na. (16) 


Letting @ = dp_i/n, Qn_1 = 0. With this value, (d,1/n) for a, a 
second special case is obtained where all the biquadratic sections (in- 
cluding the first) will have a pole-Q value equal to wo/B-n/(dn_1). For 
simple symmetrical bandpass filters, this special case reduces to 
Hurtig’s prs configuration, and Hurtig’s design formula? is identical 
to that given by eq. (3b). 

Since an infinite number of realizations, depending upon the choice 
of a, can be obtained for the shifted-companion-form representation, 
is there an optimal choice of a? This optimal choice would, perhaps, 
depend also upon the performance criterion chosen. Laker and Ghausi 
have proposed an optimization scheme for their configuration based 
on a minimization of a certain statistical multiparameter sensitivity 
measure.*:4 Their scheme can also be used here for the determination 
of an optimal a with respect to their performance criterion. In the 
following, we present two observations based on our limited design 
experience with bandpass filters using the proposed shifted-companion 
form where minimizing the filter’s passband sensitivity is of primary 
concern. In our discussion here, the filter designs are subjected to a 
computer-simulated Monte-Carlo analysis and sensitivity is examined 
from the standpoint of standard deviation (dB) vs frequency. 
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Fig. 3—Biquadratic sections Q vs shift parameter a. 


(c) It appears that a broad range of values exists for a where the 
improvements* over the cascade biquad design are relatively 
constant. This range includes Hurtig’s design, i.e., a=dn_1/n. 

(iz) Performance of the standard companion-form (i.e., a = 0) de- 
sign is about the same as that of the cascade biquad design. 


The above empirical rule (7) is observed in the design examples to 
follow. 


IV. DESIGN EXAMPLES 


Two examples given here illustrate the shifted-companion-form de- 
sign technique as well as demonstrate its performance relative to that 
of the cascade biquad, coupled (leapfrog) biquad as well as to passive 
ladder designs. Comparisons among these designs are based on a 
Monte-Carlo analysis of the filters with passive components selected 
randomly from a uniform distribution within a given tolerance interval. 


“Improvement is to be broadly interpreted as less sensitive or having a smaller 
standard deviation. 
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4.1 Example 1—-A three-section Butterworth bandpass filter* 
The normalized transfer function of a third-order low-pass Butter- 
worth filter is given by 
Vout ( ) de 1 J 
| Vin PP Op pel 
Let the desired bandpass filter have center frequency (fo) of 1 Hz 


and 3-dB bandwidth (B/27) of 0.04 Hz. The prs version of the shifted- 
companion form is designed here. Hence, 


dn-1 2 
— 7 3 





(17) 





From eq. (3a), we obtain 


a. = 0, a1 = 0.66666667, Qo = 0.25925926 
b. = bi = 0, bo = 1. 


And from eq. (5a), 
Bi = 1, B. = 0, B63; = 0. 


For this simple bandpass filter, the output summing amplifier (Fig. 2) 
is not needed. Furthermore, 


() = 008s gk 
aT ;(s) ~~ 32 +°0.087(2)s + (27)? = 1, 2,3. 
Note that 
O22 237s) Hor: 4S 4, 0:8 
*™ 0.087 2 on 


For this example, each of the 7';(s) is chosen to be realized by the single- 
amplifier biquad (sas) configuration of Ref. 8. The complete con- 
figuration’ is shown in Fig. 4, with the element values tabulated in 
Appendix B. The element values as well as circuit topologies for the 
cascade SAB, coupled sas (or leapfrog saB),* and the optimized Laker- 
Ghausi design’ are also given in Appendix B. Each of the three bi- 
quadratic bandpass sections in the shifted-companion-form, Laker- 
Ghausi, and coupled-biquad designs has a pole frequency of 1 Hz; 
whereas for the cascade design, the pole frequencies are 1, 1.01747, 
and 0.982828 Hz. The pole-Q values for these four designs are tabu- 


*'This example can also be found in Ref. 3. 

The inverting amplifier A2 can be eliminated by feeding the output of section 3 to 
the positive input terminal of the summing amplifier A:. This has not been done in 
the example. 

For symmetrical bandpass filters derived from an all-pole low-pass prototype, the 
coupled biquad and the leapfrog designs can be made identical. 
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SAB SAB SAB 
Vin(s) SECTION 1 SECTION 2 SECTION 3 





Vout(s) 


WHERE 


SAB 
SECTION i 





Fig. 4—Configuration for the three-section Butterworth bandpass filter. 


lated in Table I. These four realizations of the Butterworth filter as 
well as the passive ladder realization were compared by a Monte-Carlo 
study (with 200 trials) using the computer program BELTaAP.® The 
following assumptions are made: | 


(t) The operational amplifiers are ideal. 
(22) All passive components have the same tolerance with a uniform 


distribution. 
Table 1— Pole-Q values for the three-section Butterworth filter 
Section 
ce Pole-Q 

ilte e 

nan 1 2 3 
Shifted-companion form (PRB) 37.5 37.5 37.5 
Laker-Ghausi (FLF) 44.2 44,2 28.8 
Cascade biquad 25.0 50.0075 50.0075 
Coupled biquad/leapfrog 25.0 00 25.0 
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Fig. 5—Simulated variations of the three-section Butterworth bandpass filter 
(0.1732-percent passive components tolerance). 


Two different tolerances were simulated, the first having a realistic 
tolerance of +0.1732 percent and the second a large tolerance of 
+1.732 percent.* The resulting comparisons based on the standard 


* Realistic in the sense that the statistical variation of the filter response is within 
a reasonable bound from the nominal. The large tolerance corresponds to a component 
standard deviation of 1 percent, which was used in the example of Ref. 3. 
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deviations of the transfer function (dB) of the various designs plotted 
vs frequency (Hz) are shown in Figs. 5 and 6. 

It is observed that, over most of the passband (between the 3-dB 
points), the coupled biquad/leapfrog, Hurtig’s prs/shifted-companion 
form, and Laker-Ghausi FLF designs show roughly the same improve- 
ment (3-to-1 reduction in standard deviations) over the cascade 
biquad design. The passive filter is, however, seen to be the least 
sensitive. 


- I 
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Fig. 6—Simulated variations of the three-section Butterworth bandpass filter 
(1.732-percent passive components tolerance). 
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4.2 Example 2—A three-section elliptic bandpass filter* 


A sixth-order elliptic bandpass filter with center frequency (fo) at 
2805 Hz, 0.1-dB passband ripple with bandwidth (B/2x) of 90 Hz, 
and a minimum 30-dB loss below 2694.8 Hz and above 2919.8 Hz is 
desired. The corresponding third-order normalized low-pass prototype 
transfer function is given by 


Vevat ( ) = N(p) 
Vin D(p)’ 





where 
N(p) = 0.214115(p? + 8.158500) 
D(p) = p? + 1.897376p? + 2.543168p + 1.746858. 


Once again, the Hurtig criterion is used for the shifted-companion-form 
design. Hence, 


Q. 


2 = 0.6324587. 


a= 


oo| 


From eq. (8a), 


ad, = 0, a; = 1.8431556, ao = 0.64438116 
be = 0.214115, by = — 0.2708378, bo = 1.8325041. 


And from eq. (5a), 
B1 = 1.5449143, B2 = — 0.2708378, B; = 0.214115. 


The feed-forward zero-forming configuration (Fig. 2a) is used for the 
realizationt where 


1807s 


Pils) = 3 180n0.6324587)s + (n-2805)? 





For this example, each of the 7';(s) is chosen to be realized by the 
three-amplifier biquad configuration. The complete configuration? is 
shown in Fig. 7 with the element values tabulated in Appendix C. 
The element values for the cascade biquad and coupled biquad designs 
are also given in Appendix C.5 A leapfrog design is also available," the 
performance of which is similar to but slightly inferior to the coupled 
biquad design. Once again, a Monte-Carlo study was made on these 


* This example can also be found in Ref. 6. 

tIt was found, for this example, that the design with the feed-forward zero- 
forming technique outperforms the summation zero-forming technique design. 

+ With the three-amplifier biquad sections, it is possible to eliminate the input 
summing amplifier A; by using node 1 of section 1 as the summing point. This has 
not been done in the example. 

§ Without the availability of a computer program to choose an optimized Laker- 
Ghausi circuit, no FLF design is included in this example. 
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Fig. 7—Configuration for the three-section elliptic bandpass filter. 


designs as well as a passive ladder design. The resulting comparisons 
are shown in Figs. 8 and 9, where a 0.25-percent tolerance (with 
uniform distribution) is assumed for all passive components. 

It is observed that, in terms of standard deviation, the passband 
improvements over the cascade design are noticeably less for the 
shifted-companion-form (a@ chosen by Hurtig’s PBR criterion) design 
than the roughly 4-to-1 improvement of the coupled biquad design. 
Once again, the passive filter outperforms its active counterparts. 


V. CONCLUSIONS 


The FLF/PRB multiple-loop feedback active filter structure is known 
to have better sensitivity performance than the popular cascade ap- 
proach. This sensitivity improvement is particularly acute in high-Q 
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Fig. 8—Simulated variations of the elliptic bandpass filter. 


bandpass filter designs, as exemplified in the two examples given. With 
the described shifted-companion-form representation of the filter 
transfer function, it is straightforward to obtain explicit design formulas 
for this feedback structure as contrasted to the coefficient matching 
technique used by Laker and Ghausi. In practice, the shift parameter 
can be chosen such that identical biquadratic blocks (i.e., the extended 
PRB version)* are used. The proposed shifted-companion-form design 
has the following advantages over the optimized Fir design. First, 
no matrix inversion and involved sensitivity minimization routine 
are needed. Second, all sections are identical, and the sections’ pole-Q 
can be much lower than the highest pole-Q required for the FF design. 
Furthermore, little difference is usually observed for this shifted- 


* Hurtig’s pRB does not treat the cases with finite transmission zeros. 
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Fig. 9—Simulated variations (passband) of the elliptic bandpass filters. 


companion-form design and the optimized Fur design. On the other 
hand, the two examples also show that the coupled biquad and leap- 
frog designs may have better sensitivity performance than the cor- 
responding shifted-companion-form designs. However, for those types 
of filter functions having finite transmission zeros, the designs of 
coupled biquad and/or leapfrog require extensive computer aids 
that are not yet generally available. Hence, the proposed shifted- 
companion-form technique may provide an alternative to the coupled 
biquad/leapfrog active filter design techniques. 
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APPENDIX A 
Derivation of the Shifted Transfer Function 


Let the polynomial P(p) be 


P(p) = Vides (18) 
and let 

p=s—-a. (19) 

Then 

P(s) = > dn—i(S — a"? 

7 = dn: & . :) seria), (20) 

where 

ra. (aia)! 

( r ) m rin—i—r)l (21) 


Equation (20) can be rearranged in decreasing power of s as follows: 


P(s) = d,s* + 8" 3 « 7 ") (—a)'~"d.| 


r=0 


icenens 
$e port (ETT) (aed 


r=0 \2 


ie n—r — 
+: +E (N TT) He) ae 


Or 
P(s) = > gn-t > (—1) ( es ) adn 
i=0 r=0 1-?r 
= 3 an_ss" 7, (22) 

i=0 

where . 
: wr ( OT YS eee 
ani = PDA Get € = ya is, 

or 


> a irq. 
r=0 ( 


iG (23) 
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Equation (23) can be used to obtain the coefficients of the shifted 
polynomial P(s) from the coefficients of the polynomial P(p). 


Similarly, if we start with the polynomial 


P(s) = 3 Gast 


7=0 
and let 
s=p+a, 
then we obtain 


P(p) = ¥ dap", 


where , 
: n—r)! so 
c= 2G pGa ale Oe 


APPENDIX B 
Element Values for Example 1 (Fig. 4) 


(24) 


(25) 


(26) 


For the various realizations, resistors are in kilohms and all ca- 


pacitor values are 10 pF. 


B.1 Shifted-companion-form (Hurtig’s) realization 


Section 
ieee 1 2 

R2 128.5 128.5 
Ry 613.2 613.2 
Rs 1.977 1.977 
Ry 73.07 73.07 
Ra 2.0 2.0 
R, © re) 


and Rin = 2.963, Rye = 6.667, Ry = 11.48. 


B.2. Laker-Ghausi realization 


Section 

Blemen< 1 2 
Re 128.5 128.5 
Ra 722.3 722.3 
Rs 1.976 1.976 
Ry 71.77 71.78 
Ra 2.0 2.0 
R, ee) °o 


and Rin = 2.781, Ryz = 4.600, Rys = 23.73. 


SHIFTED-COMPANION 


128.5 
613.2 
1.977 
73.07 
2.0 


127.7 
470.9 
1.993 
74.73 
2.0 
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B.3 Cascade realization 


The summing amplifier A; and inverter Az are not needed, and the 
input goes directly to node 1 of section 1. 


Section 

Sona 1 
R, 128.5 
R, 408.1 
Rs 1.981 
R, 77.76 
Ra 2.0 
R, ro) 


B.4 Coupled biquad realization 


126.3 
409.1 
1.946 
70.93 
2.0 


130.8 
409.1 
2.015 
70.93 
2.0 


The summing amplifier Ai and inverter A» are not needed, and the 
input goes directly to node 1 of section 1. 


Section 1 
Element 


R. 128.5 
R, 204.1 
Rs 1.990 
Ry 77.39 
Ra 2.0 
R, 414.1 


2 


128.4 
820.2 
1.977 
64.98 
2.0 
829.8 


128.5 
408.1 
1.981 
77.76 
2.0 


In addition, node 1’ of section 1 is connected to node 2 of section 2 and 
node 1’ of section 2 is connected to node 2 of section 3. 


APPENDIX C 


Element Values for Example 2 (Fig. 7) 


For the various realizations, resistors are in kilohms and all ca- 


pacitor values are 0.01 uF. 


C.1 Shifted-companion-form (Hurtig’s criterion) realization 


Ri 279.6 
Ro 5.674 
R; 5.674 
Ry 279.6 
R; eo 
R, eo 

Ry 10.0 
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2 


279.6 
5.674 
5.674 

279.6 

eo 


© 
10.0 


279.6 
5.674 
5.674 

279.6 


eo 
10.0 


1975 


and 
Rye = 2.978, Ry = 3.926 
Rint = 1.638, Rine = 413.0, Ring —) ae 825.9.* 


C.2 Cascade realization 


The summing amplifier A; and all input feed-forward paths (Ri,’s) 
are not needed. The input goes directly to node 1’ of section 1. 


Section 
aie 1 2 3 

Ry 167.5 428.6 412.2 
Ro 5.674 5.786 5.564 
R; 5.674 5.786 5.564 
Rs 167.5 1190.0 732.2 
Rs 00 16.93 9.378 
Re re) 27.75 17.76 
Ry 10.0 10.0 10.0 


C.3 Coupled biquad realization 


The summing amplifier A, is not needed. The input, Vin, goes 
directly into node 1’ of section 1 and nodes 1 of sections 2 and 3 
through the feed-forward resistors Ring and Rins, respectively. 


Section 
Element J 7 3 


Ri 311.5 00 133.0 
Re 5.674 5.674 5.674 
Rs 5.674 5.674 5.674 
Ie4 97.92 oO re) 

Rs eo oe} oo} 
Rs eo lee} oe} 
Ry 10.0 10.0 10.0 


and Ring = 1324.0, Ring = 825.9. In addition, the following resistors 
are needed with value and connections noted. 


(t) 180.5 kQ, node 2 of section 1 and node 1 of section 2. 
(zz) 180.5 kQ, node 1 of section 1 and node 4 of section 2. 
(277) 194.3 kQ, node 2 of section 2 and node 1 of section 3. 
(tv) 194.38 kQ, node 1 of section 2 and node 4 of section 3. 


*In practice, with the following modifications of Fig. 7, Rin3 = 825.9 is used. 
Change the connection of sections 2 and 3 to between node 2 (section 2) and node 1’ 
(section 3); the connection of Ryz remains unchanged. Change the connection of 
Ry3 to between node 2 (section 3) and the summing amplifier A. 
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Timing Recovery and Scramblers 
in Data Transmission 


By R. D. GITLIN and J. F. HAYES 
(Manuscript received August 28, 1974) 


This paper considers several problems associated with envelope-derived 
timing recovery, equalization, and scrambling in synchronous data trans- 
mission. Particular attention is focused on the time intervals in which 
periodic data sequences are transmitted, such as during start-up or when 
an idle code is being transmitted. It is shown that the standard envelope- 
derived timing-recovery system may be significantly improved by zonal 
jiltering of the received passband signal prior to forming the envelope. For 
phase-modulated systems, we discuss the limitations of the “precession” 
technique employed for the purpose of providing a pertodic timing wave 
when there is an input of short period. The advantages of using a phase- 
locked loop to filter the envelope instead of a narrow-band filter are also 
described. A study of scrambler operation has provided an extension of 
previous results concerning the relationship between the input and output 
period. It 1s shown that the output period of several scramblers connected 
an tandem does not necessarily double with the addition of a stage, and 
that af a particular stage does not lock up then no succeeding stage can 
lock up. 


I. INTRODUCTION 


Recovery and tracking of the symbol rate, or timing frequency, is 
one of the most critical functions performed by a synchronous modem. 
Most modems are “self-timed” in that they derive their timing fre- 
quency and phase directly from the information-bearing signal, instead 
of using a separate subchannel to send synchronization information. A 
technique that is commonly used to acquire the symbol rate* (which 
is the receiver’s basic sampling rate) is to filter the envelope of the 
modulated data signal. Our investigation will consider several problems 
related to this method of timing recovery which arise in high-speed 
modems incorporating both an adaptive equalizer and a scrambler. 


*This technique is also used to provide the sampling epoch, or phase, within a 
symbol interval. 
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The envelope-derived timing recovery system is a well-studied 
topic.1:? However, as the degree of excess bandwidth decreases, the ease 
with which timing can be recovered using this approach rapidly di- 
minishes. We focus our attention on periodic input sequences. These 
sequences are used to train (or adapt) the data receiver during start-up 
and during the idle period between blocks of random data. To provide 
a densely spaced line spectrum of uniform amplitude (which is neces- 
sary if the equalizer coefficients are to remain properly adjusted for 
random input data), high-speed modems use a scrambler to ‘‘random- 
ize’”’ the short periodic inputs commonly used during the idle period. 
We investigate the effect of the scrambler on both the line spectrum 
and the strength of the timing tone. It is observed that zonal filtering 
of the received data signal prior to taking the envelope can signifi- 
cantly improve the relative strength of the timing tone by suppressing 
the jitter component. 

Using transform theory, a discussion is presented on the relation- 
ship between the scrambler input and output periods. We refine 
Savage’s? well-known results for periodic inputs; these refinements 
are applied to the study of the tandem and parallel scrambler 
configurations. 

Sections II to IV review the basic envelope-derived timing system 
and give expressions for the power in the timing and interfering tones. 
The role of the phase-locked loop in the timing recovery system is 
described in Section V. Section VI considers the effect of precessing” 
the data symbols on timing recovery. The necessary background 
material on self-synchronizing scramblers is presented in Section VII. 
The transform approach is used in Section VIII to determine the 
scrambler output period. In Section IX the performance of a cascaded 
scrambler configuration is contrasted with the conventional serial 
arrangement. The parallel scrambler configuration is discussed in 
Section X. 


Il. BASIC TIMING RECOVERY SYSTEM 


In this section we describe the commonly used technique of acquiring 
the timing frequency and phase by processing the envelope of the re- 
ceived signal. The object is to extract a tone, located at the symbol 
rate, which is then used in the sampled-data receiver. Figure 1 shows 
a simplified receiver structure of an in-phase and quadrature (e.g., 
QAM) data-transmission system, where we have focused attention on 
the timing recovery and equalization functions of the receiver. For 


* The advancing of the transmitted angle by a fixed phase (in a differential phase- 
modulated modem), independently of the input, is known as precession. 
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our purposes it will be convenient to ignore both the additive noise 
and the quadrature component of channel distortion. Using the nota- 
tion of Figs. 1 and 2, the received signal s(t) is given by 


{3 } DETECTED 

2 DATA 
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meee 

ENVELOPE PHASE-LOCKED 
LooP 


Fig. 1—Simplified Qa receiver. 











st) = D> ang — nT) coswt — YO dbag(t — nT) sinw.t, (1) 
where a, and b, are respectively the discrete-valued in-phase and 
quadrature data sequences obtained from the binary sequences an 
and Bn, g(t) is the spectral-shaping pulse, w, is the carrier frequency, and 
1/T is the symbol rate or timing frequency. The envelope of a filtered 
version of the received signal is tracked by a phase-locked loop tuned 
to the receiver’s best a priori knowledge of the timing frequency. The 
output of a properly designed phase-locked loop will be a tone with 
frequency equal to the symbol rate and whose zero crossings may be 
used to derive a sampling wave. Once the timing frequency is acquired, 
the estimated and unscrambled data sequences {@,} and {8,} are 
available to the user. The decoder maps the sequence of multilevel 
two-tuples (dn, 6») into a binary sequence which serves as the input to 
the inverse scrambler. 

The data sequences {an} and {b,} may be thought of as random 
(when user data are being sent) or as periodic (during start-up when 
the equalizer and timing parameters are being acquired, and during an 


COS wet 
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INPUT TRANS— 
BINARY MITTED 
SEQUENCE SIGNAL 


FILTER 


FILTER 


SIN wet 


Fig. 2—Simplified Qam transmitter employing a scrambler. 


TIMING RECOVERY AND SCRAMBLERS _ 571 


idle period between random data transmissions). As we shall shortly 
see, the presence of a short periodic input can play havoc with the 
equalizer tap settings; hence, a scrambler is generally used at the 
transmitter to ‘‘randomize”’ these periodic sequences. As depicted in 
Fig. 2, the coder maps the binary stream of scrambled Os and 1s into 
the channel pulse levels (e.g., 0 and 1 may be mapped into —1 and 1 
respectively).* The effects of choosing a particular scrambler structure 
(e.g., serial vs parallel or serial vs cascade) on the timing recovery 
system will be treated in Sections VII and VIII. 


ill. SPECTRUM OF THE RECEIVED SIGNAL 


We confine our attention to periodic inputs, beginning with a cal- 
culation of the Fourier transform of the received signal. Rewriting (1) 
in complex notation, we have 


s(t) =Re} SY cag(t — Teel , (2) 


where cn = Qn + jb, and Re denotes the real part of a complex number. 
With a periodic data sequence, cn, the signal s(t) is periodic. This latter 
periodicity is best exhibited via the discrete Fourier transform! (pFr) 
of the periodic sequence. With the period of c, denoted by N, the prr 
of c, is defined by 


N—-1 
CRD): = SY egen teat k=0,1,:°-,N-1 (8a) 
n=0 
and the inverse relation is 
1 N=! 
Cn = = DY Clk) esn0r n=0,1,---,N —1, (3b) 
N £=0 


where 2 = (1/N)(2x/T) = (1/N)- (symbol frequency). Hence, the 
prt has N components uniformly spaced 1/NT Hz apart and the spec- 
trum repeats every 27/7’ Hz. Denoting the Iourier transform of s(¢) 
by S(w) and convolution by @ we have 


S(w) = 4 LD cre" TG(w) |®S(w — «-), w > 0, (4) 
and using (3b) in (4) givest 


Lt 2rn 27n ; 

Sto) = 5 2, Cem) | Ea (na + °F) 5(o— a. — na — 7) 
w>O0. (5) 
Letting the timing frequency be denoted by w, = 27/T = NQ, it is 


*In practice, the data would also be differentially and Gray encoded. 
t We use the identity >, e"*? = (27/T) on d[w — (2an/T)]. 
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W- > 
Fig. 3—The spectrum of the line signal for a periodic input is 
S(w)=Re {81 C(kQ) [4 (kQ)5 (w —we —k2) + E(k —ws)8(w —we—kQ+ws) ]} 
where the input period is NT’ seconds. 


clear that S(w) has discrete tones at w,. + kQ + nus. In practical data- 
transmission systems with pulses of less than 100-percent excess band- 
width, G(kQ + 21n/T) will be zero if n ¥ 0 or —1, hence, 


SW) = 57 ps C(kQ)[G(k2)5 (eo — w. — bQ) 
+ G(kQ — ws)6(w — wo, —kQ—s)]; w>O0. (6) 


This spectrum is illustrated in Fig. 3, where it is seen that the envelope 
C(kQ) modulates the baseband pulse shape, G(w — w-), in the range 
We — We tO We + ws. 

Since the signal s(¢) is used by the equalizer to adjust the tap weights, 
ideally the spectrum C(kQ) should approximate that of random data, 
i.e., be constant. Of course, it is more critical that the equalizer be 
presented with a closely spaced line spectrum; for example, if the 
input period were two symbol intervals, it is clear that the equalizer 
can only compensate for distortion at two frequencies in the Nyquist 
band. Consequently, at the instant when the data return from the 
periodic to the random mode, the equalizer tap settings will be far from 
their optimum (for random data) values, and the distortion at the 
equalizer output could be much larger than the channel distortion. This 
situation generally causes the receiver to make so many errors that it 
is necessary to retrain the equalizer. As we shall see, the role of the 
scrambler is to lengthen the period of the transmitted sequence, 
thereby keeping the equalizer trained. Hence, for the rest of our dis- 
cussion, we will assume that the scrambler is such that the periodic 
spectrum is (essentially) flat and densely spaced. Section VIII deals 
specifically with the factors that determine the period of the scrambler 
output. 


IV. SPECTRUM OF THE ENVELOPE 


The timing frequency is to be acquired from the envelope of the 
filtered line signal. Let the filtered line signal m(é) be 


mt) = Yo anf(t — nT) cosad¢ — d baf(t — nT) sin wd, (7) 
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where f(é) is the (equivalent baseband) pulse shape after filtering at the 
receiver. The (squared) envelope of m(t) is defined as 


rt) = LX anf(t — nf) P + 0X daft — nP) P. (8) 


As before, we introduce complex-signal notation by letting 


d(t) =  enf(t — nT) 
BO) = Last — v7), (9) 


so that we can write 
r(t) = d(t)-d*(@), (10) 


where * stands for the complex conjugate. Thus, the Fourier transform 
of r(t) is given by 


R(w) = D(w)@D,(@) = [Li cre" (wo) JOE Cre"? F(w)], (11) 


where D,,(w) is the Fourier transform of d*(¢), and F'(w) is the transform 
of f(t). Using (8b), we have that 


D(w) 


N-1 
y | Xs Cae serr>¥ e— sent F (wy) 


N-1 

~ C(kQ)F(w) ¥ 6(w — kQ — nas). (12) 
k=0 n 

Substituting (12) into (11) and performing the convolution gives 


ene = ba C(k9)C*(10) © F(KQ + nox) F (mos — 19) 


X s[w — (k —DQ— (n+ mo]. (18) 


Evidently there are tones at pQ+qw; (where p=k—I and g=n+m); 
the desired tone is at ws (i.e., p = 0, g = 1) and all other tones may 
be regarded as interferers. Again, practical bandlimiting of F'(w) and 
filtering of R(w) will eliminate all terms where q ~ 0 or 1. The power in 
the desired tone is 


N—-l 
R(ws) = D> |C(kQ) |PF(kKQ)F (ew, — kQ), (14) 
k=0 
while the power in an interfering tone (or sidelobe) is 


R(w. + p9) = 3 C(kQ)C*E(k — p)9]F(b9) FC (k — p)@ — ws] 
p= Ty 2, 3S" (15) 
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Fig. 4a—Strength of timing tone without loop filtering is 
UIC (KQ) |?G(kQ)G (ws — KO). 


G ( Ws—(k—p).N), p>0 






u-a)5s = (1 +a 
Fig. 4b—Strength of sidelobe interference at w, + pQ is 
Litt C(kQ)C*L(k — p)QG(kQ)GL(k — p)Q — we]. 


Conventionally, the signal r(¢) is fed to a phase-locked loop (PLL), 
which acts like a narrow-band filter in accepting tones in some region 
about ws (e.g., from w, — BQ to w. + BQ, where 2BQ is the effective 
bandwidth of the loop) and produces an output that is dominated by 
a tone at ws. We first consider the above spectra in the absence of any 
timing loop filtering [i.e., F(w) = G(w) ] as shown in Fig. 4. It is clear 
from (14) that the problem of timing-frequency recovery becomes more 
difficult as the amount of excess bandwidth (as measured by the 
parameter a) decreases—for zero excess bandwidth, this timing re- 
covery technique clearly fails since the pulses G(kQ) and G(w, — kQ) 
are disjoint. Figures 4a and 4b show how to compute the power of the 
tones at w, and at ws + pQ respectively. We note that, in general, 
R(w, + pQ) # R(w. — pQ), and moreover, for the particular spectral 
shaping shown in the figure, it is clear that R(w. + p12) > R(ws + po) 
for —B < po < pi < B; thus, half of the sidelobe tones are greater in 
magnitude than the desired tone. Thus, without any prefiltering in the 
timing loop, the desired tone is rather weak in comparison to the 
interfering tones. As we have already mentioned, this problem has a 
direct solution’: choose the loop filter F(w — w.) to be a narrow zonal 
filter around w, + (w;/2) and w, = (w;/2) as shown in Fig. 5. The 
resulting signal m(¢) has its energy concentrated at w, — (w,/2) and 
We + (ws/2), and the relative strength of the timing tone is illustrated 


“A filter in the timing loop has also been proposed by Franks and Bubrowski.® 
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F(w—we) 


@ 
Ws We Ws 
Wc 2 Wet 2 


Fig. 5—Timing-loop filter shape which improves the tone-to-interference ratio. 


in Figs. 6a, 6b, and 6c. The attenuation of the interferers is aided 
further by the fact that the magnitude of -¥=j C,.Cy_, is a maximum 
for p = 0 (this follows from the Schwarz inequality). Clearly, by 
making F(kQ) = 6(kQ — ws), we can make R(w, + pQ) = 0 for all 
p ~ 0; however, any narrow-zonal prefilter of the type shown in Fig. 5 
should significantly improve the relative strength of the timing tone. 
Since we merely require the filter to be narrow-band, any reasonable 


ka 





Fig. 6b—Strength of tone at w, + pQ ~ LNG CiCy_,F (kQ)F E(k — p)Q — we]. 


ka 
Ws 


Fig. 6c—Spectrum of envelope. 
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LOW-PASS e(t} 
FILTER 






COS(wst +f e(t)dt) 
Fig. 7—First-order phase-locked loop. 


choice wide enough to accommodate the uncertainty in w, would be 
adequate over a wide range of channel characteristics. 


V. THE ROLE OF THE PHASE-LOCKED LOOP 


As we have discussed, the signal r(¢) contains a desired tone as well 
as interfering tones. In this section, we wish to demonstrate that a 
phase-locked loop (PLL) can provide extremely narrow-band filtering 
even to the extent of extracting a single tone from a spectrum of 
adjacent interferers. Consider the standard first-order PLL® shown in 
Fig. 7. Let us assume that the input is the desired tone plus two 
interfering sidelobes, i.e., 


r(t) = Asinwst + Bsin[(w, + At+ y] 

+ Bsin[(w — A)jt-—y], (16) 
where A is the frequency displacement of the sidelobe from the desired 
tone and y is the corresponding phase shift. Note that we have special- 
ized the situation to the case where both interferers have the same 
amplitude and opposite phase angles (i.e., the distortion in the timing 
recovery system is symmetric about w, radians). We also assume that 
a perfectly tuned loop (i.e., the free-running frequency of the voltage- 
controlled oscillator (vco) is w,) is employed.* From Fig. 7 the loop 
error signal is given by 


e(t) = Asin ( f e(at + a) + Bsin( fe(dt +a — At— 7) 
+ Bsin( [edt +a+ att 7): (17) 


If we define ¢(¢) as the phase difference between vco output phase and 
the PLL input phase corresponding to the desired tone, i.e., 


o(t) A wt — (os + / e(t)dt + a), (18) 


“A perfectly tuned loop could arise by varying the free-running vco frequency. 
As we show, via (22), when this condition is achieved the output will be a single tone. 
This observation suggests a feedback or error-sensing procedure for varying the 
nominal vco frequency. 
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it is necessary that ¢(¢) 0, as t+, since this implies successful 
tracking of the tone. Using (18) in (17) gives 


e(t) = — [A + 2B cos (At + y) ] sin o(d), (19) 
and since from (18) ¢(t) = — e(t), the PLL is governed by the first- 
order differential equation 

ew) 2 AOR aw Cnn ia tease: (20) 
To solve (20), we first separate the variables and write 

dp _ | 
ne ae A+ Bcos (At + y)dt, (21) 


and, by direct integration, we obtain the solution 
o(t) = 2tan— {e-4* exp (2B/A)[sin (At + 7) — siny]}. (22) 


We then have ¢(t) ~0 as t->~, i.e., the loop locks on the desired 
tone for any strength of the interference tone. Clearly, the same would 
be true for a collection of interferers provided they met the assumed 
symmetry conditions on their amplitude and phase. This example 
illustrates the power of a PLL to capture a desired tone in the presence 
of considerable interference. 


vi. EFFECT OF PRECESSION ON THE RECOVERY OF A TIMING TONE 


In modems not employing adaptive equalization, the question arises 
as to whether or not a scrambler is needed to generate a timing tone 
during the idle period. Since there is no equalizer in the system, we 
are not concerned with having a densely spaced line spectrum but only 
that there be at least two spectral lines, spaced w, apart, in the pass- 
band signal. Using the framework we have developed in the preceding 
sections, we investigate the effect of ‘“‘precessing”’ the data symbols. 
Let us consider the phase-modulated signal 


s(t) = ps g(t — nT) cos (wet + On), (23) 


whose idle code is 6, = 0 for all n. The spectrum of s(¢) is, by using 
(5) with C(kQ) = OKO and N => 1, 


S(w) = = CO io. (24) 


and for an excess bandwidth of less than 100 percent, 


S(w) = G(0)d(w — a) ; (25) 
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fo = 1800 Hz 
1/T = 1200 SYMBOLS / SECOND 
100% EXCESS BANDWIDTH 


f (Hz) 





1200 1800 2400 3000 


Fig. 8a—Spectrum without precession. 


f (Hz) 
750 1950 


Fig. 8b—Spectrum with precession. 


i.e., the spectrum consists of a single tone, which is obviously not 
sufficient to provide a timing tone. This situation is illustrated in Fig. 
8a for a pulse g(t) with 100-percent excess bandwidth and with f, = 1800 
Hz and 1/T = 1200. To avoid the above situation, let 6, = 2n7/M, 
where J/2 is the number of points in the signal constellation and, thus, 
6, has period M. The advancing of the data symbol by 27/M degrees, 
independently of any change in the input data, is known as precession.* 
Using the notation in Section I, we have 


Cn — ei6n = es nr] M) 


M-1 : MeN.) 2. 
CEQ) = 2 cnc iabarsa = YY g-i@risnat—) = 5, (26) 
n=0 


n=0 


which from (5) gives 
S(w) = LG (nw, + Ff) 3(o — we — GF — nw.) (27) 


Thus, the effect of precession is to offset the tones by w,/M, producing 
the spectrum shown in Fig. 8b. Clearly, when squared, this signal 
provides a tone at w, = 1200. 

The situation is different, however, for the spectrum shown in Fig. 9a 
with 6, = 0. The spectrum with precession shifts the tone by 100 Hz, 
and clearly no pair of in-band tones is present. Thus, for spectra that 


* Differential phase modulation with precession would generate a data sequence 
On = On-1 + bn + 2n7/M, where $n is one of M/2 equally spaced angles. 
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50% EXCESS BANDWIDTH 
fo = 1800 Hz 
1/T = 1600 


f (Hz) 


200 600 1800 3000 3400 


Fig. 9a—Signal spectrum without precession. 


f (Hz) 
300 600 1900 3000 3500 


Fig. 9b—Signal spectrum with precession. 


use small amounts of excess bandwidth, the tones necessary to provide 
timing would not be present with or without precession. Also, preces- 
sion has little or no bearing on equalizer training since it simply shifts 
the spectrum and does not enrich it. For high-speed transmission in 
which the amount of excess bandwidth is small, spectral enrichment is 
provided by the scrambler and insures the proper operation of the 
equalizer and timing recovery system. 


Vil. SCRAMBLERS: BACKGROUND MATERIAL 


We have shown in the preceding sections how the transmission of 
short repetitive patterns can play havoc with both the equalizer and 
timing recovery systems. As the name suggests, scramblers serve to 
“randomize” deterministic data sequences. The effect of this random- 
ization on periodic sequences is to lengthen the period of the input 
sequence to the scrambler. Strictly speaking, the periodic output of 
the scrambler is not random. However, the scrambled data stream 
results in a line signal that has many more spectral components than 
the input data stream, and, thus, it looks more like the continuous 
spectrum that results when purely random data are encoded. 

We confine our attention to the so-called self-synchronizing scram- 
bler.2 The generic forms for the self-synchronizing scrambler and the 
descrambler are shown in Figs. 10 and 11 and consist of, respectively, 
feedback and feedforward shift registers. Data symbols are fed into 
the scrambler every T seconds. These symbols are added (modulo p)* 


“In practice, p = 2. 
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m 
Yn = Un © bs) Ci Vn-i 
i=1 








Fig. 11—Inverse scrambler. 


to past outputs to produce the current output. The inputs to the delay 
elements shown in Fig. 10 are delayed by 7 seconds. The output of the 
scrambler is then encoded for transmission over the channel. After 
decoding at the receiver, the resulting sequence is put through a de- 
scrambler, shown in Fig. 11, where the original sequence is recovered. 
The inverse scrambler is self-synchronizing, and it will eventually 
cleanse itself of a transmission error once the error has propagated 
through the shift register. The number of errors in the descrambler 
output sequence is the number of channel errors multiplied by the 
number of nonzero tap weights in the shift register. 

We shall study the input-output relationships of scramblers using 
d-transforms. Using this tool we are able, quite simply, to extend and 
clarify Savage’s theorem’ on scramblers with periodic inputs. Before 
getting into details on scrambler input-output relationships, a sum- 
mary of some necessary background material on polynomials over 
Galois fields is in order.* 

With p a prime number, we speak of a polynomial Q(d) over GF (p), 
where the coefficients of Q(d) are modulo-p numbers. Multiplication, 
addition, and division of such polynomials are carried out in the usual 
fashion using modulo-p arithmetic on the coefficients. The degree of a 
polynomial Q(d) is the highest power of d appearing in Q(d). A poly- 
nomial of degree m is irreducible if it cannot be factored into poly- 


* Much of the background material is taken from Ref. 7. 
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nomials of lower order. Two polynomials are relatively prime if they 
have no common factors. A crucial concept in our study of the scram- 
bler is the exponent of a polynomial. The exponent of the polynomial 
Q(d) is the minimum value of | such that Q(d) divides 1 — d’, ice., 
(1 — d')/Q(d) is a polynomial of finite degree. For example, the 
exponent of the polynomial 1 + d? + d* in GF (2) is 7 since it divides 
1+ d’, yielding 1+ d?+ d?+ d‘, but it does not divide 1 + di‘, 
«<7. If the polynomials P(d) and Q(d) are relatively prime with 
exponents J; and J: respectively, it can be shown that the exponent of 
P(d)Q(d) is the least common multiple (lem) of 1; and lz. The ex- 
ponent of [Q(d) ]’, where Q(d) is over GF (p), is p’l, where I is the ex- 
ponent of Q(d) and r is such that p™ <j S p*. An irreducible 
polynomial of degree m is primative or of maximum exponent if its 
exponent is p” — 1. Given a polynomial Q(d) of order m, its reciprocal 
polynomial is d”Q(1/d), and it is known that reciprocal polynomials of 
irreducible polynomials are themselves irreducible, and that reciprocal 
polynomials of primative polynomials are themselves primative. 

This theory of polynomials over a Galois field is applicable to the 
d-transforms’ of the input and output sequences of a scrambler. Con- 
sider a time series %, 21, 2, --:, such that the z,,7 = 0, 1, --- are 
elements from a Galois field, e.g., 01101---. The d-transform of this 
series is defined as 


X@ = y ad‘, (28) 


and inversion is accomplished by “reading” the coefficients of X(d). 
The d-transform of a periodic sequence is of the form R(d)/(1 — a), 
where ) is the period of the sequence and R(d) is a polynomial, of 
degree less than A, in d over GF(P). To illustrate, suppose we have a 
series of elements in GF (3), 1021, 1021, ---. Using the relationship for 
a geometric progression we find that the d-transform of this series 
is (1 + 2d? + d*)/(1 — d*). In general, it can be shown that the 
d-transform of a periodic time series is of the form P(d)/Q(d), where 
P(d) and Q(d) are polynomials over a Galois field. If P(d) and Q(d) 
are relatively prime, the period of the time series represented by 
P(d)/Q(d) is the exponent of Q(d). 

Linear sequential machines over GF(p) are composed of modulo-p 
adders, multipliers, and delay elements connected according to a few 
elementary rules. As the name implies, such circuits are linear over 
modulo-p arithmetic. The laws of commutativity, associativity, and 
superposition apply. For example, the response of a circuit to the sum 
of two inputs is the sum of the responses to each input separately. The 
summations are carried out term-by-term modulo p on the input and 
output sequences. 
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The response of such circuits to inputs can be found using the 
classical techniques of linear system theory. The output consists of the 
sum of the free response and the forced response. The free response is 
due solely to initial conditions within the circuit with no input. If the 
circuit is in the quiescent state, i.e., it has zero initial conditions, there 
is no output without an input. The forced response is the output when 
an input is applied to a circuit in the quiescent state. As in the case of 
conventional linear circuits, the forced response can be found by con- 
volving the impulse* response with the input sequence. Thus, if An, 
n = 0,1, --- is the impulse response of a circuit and wu, is its input at 
time nT, n = 0, 1, ---, then the output at time nT is 


Yr = > hiUn—r, (29) 
k=0 


where ») indicates summation modulo p. If we take the d-transform of 
both sides of (29), we find 


Y(d) = U(d)H(d), (30) 


where U(d) and H(d) are the d-transforms of ua and ha respectively, 
and where n = 0, 1, ---. 


Vill. SCRAMBLER INPUT-OUTPUT RELATIONSHIPS 


Scramblers are linear sequential circuits and their input-output 
relationships can be found using linear system theory. In this section, 
we wish to demonstrate the utility of the d-transform approach in 
characterizing the nature of the output sequence for a given input 
sequence. Consider the m-stage scrambler shown in Fig. 10 with feed- 


back coefficients ci, ---, Cm. The output at time n7, yn, 1s given by 
Yn = C18in D C282n D +++ GD CmSmn@® Un 
Sin = Yn-1 (31) 
Sin = Si-1,n-1 12 2, 


where Un is the input at time nT and s;: is the output of the kth delay 
element at time J.' Now we find the impulse response. Let 


‘pe 1 n=0 
"10 n>0 


and 
S30 = 0, Vi. 


“By an impulse we, of course, mean a time series which is unity at the time origin 
and is zero elsewhere. 

T The output may be rewritten as yn = BD ciyn-i ® un. At the descrambler we 
form Zn = ¥n ® DM, Cini, which recovers the input when there are no channel 
errors. 
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The output sequence can be written 


lA 
= 


Dicyrn1 nem 
Yo=5 (32) 


D Cini n>m. 


If we take d-transforms of both sides of (32) we find after some manipu- 
lation that the transform of the impulse response, i.e., the transfer 
function, is given by 


tS = yed? = i/(1 7 2 ca) ee) 


The d-transform of the forced response of the scrambler can be found 
from eqs. (30) and (38). 

The free response of the scrambler can also be found from (81) 
when un = 0, for all n. We begin by assuming a particular initial 
state vector. Assume that the output of all of the delay elements but 
one are zero. Let the nonzero output be that of the 7th delay element 
and denote this output by sj. It can be shown that y;,, the output of 
the scrambler due solely to state s;o, is 


Ci8i0 n=O 
n 
DW CYn—7 + Ci4nSi0 O<nim-i 
Yin — j=l (34) 
2 CiYn—7 n>m—t1. 
ie 


If we take the d-transform of both sides of eq. (84), we find that 
the d-transform of the response to initial condition so is 


Yi(d) = (si > cist) /(1 — 5p eid’). (35) 
a= j=l 
Now, to find the response to any initial condition 810, S20, +--+, Smo, 


we sum over 7. Thus, the d-transform of the free response is 
Y trea(d) = S(d) if (a =e oid’) = S(d)H(d), (36) 
j=l 
where S(d) 4 Bi, sio prt Crvxd”: 
A fact that is crucial to our analysis in the sequel is that the poly- 
nomial S(d) spans the space of polynomials of degree m — 1 in GF(p). 


By choosing the initial conditions sio, i = 1, 2, ---, m, S(d) can be any 
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polynomial of degree m — 1. To show this, suppose we have the 
polynomial T'(d) = t% +tid+ --- + tnd”. Equating T(d) and 
S(d) term by term, we have m equations in m unknowns, 810, 829, - °° 
Smo. The equations can be represented in the form 


) 


Cm 0 cack 0 S10 | 
a mi-t: O Ss bigs 

° : : “ fal : vhs (37) 
Ci C2 "t+ Cm Smo to 


The m X m lower triangular matrix in (37) is nonsingular (since 
Cm #0); therefore the m simultaneous equations have a unique 
solution. 

From the foregoing, we see that the total response of a scrambler to 
an input, with transform U(d), is 


Y(d) = [U(d) + S(d@)]/e(@), (38) 


where ®(d) 4 1 — 07, c:d™ is the transform of the feedback co- 
efficients. The above equation completely describes the behavior of 
the scrambler to any input for any given initial state. 

Now we consider the input-output relationships for the scrambler 
based on eq. (38). Throughout our discussion we shall assume that 
(d) is a primitive polynomial implying that it has exponent ¢=p”—1, 
and thus can be written as (1 — d*)/’(d), where ®’(d) is a finite 
degree (remainder) polynomial of degree ¢ — m. Note that ’(d) is 
one ‘‘cycle’’ of the periodic polynomial 1/@(d). Suppose that the input 
is zero, the transform of the output is simply S(d)/#(d). Since the 
degree of S(d) is one less than that of &(d), S(d) and ®(d) are relatively 
prime,* and the output transform is S(d)®’(d)/(1 — d*); hence, the 
output is periodic with period = p™ — 1. If the input is a sequence of 
finite duration j, then U(d) is a polynomial of degree j — 1. If j S m, 
then the above output transform is U(d)®’(d)/(1 — d*), and since 
the degree of the numerator is less than ¢, it is clear that the output 
is purely periodic with period ¢. Note that there is no output transient. 
If 7 = m, and if U(d) + S(d) and &(d) are relatively prime, then it is 
easy to show that the output consists of a transient component 
(7 + 1 — m) longt and a periodic component with period ¢. For any 
input U(d) of finite duration, there are a unique set of initial conditions 


*Since @(d) is a primative polynomial, S(d) cannot be a factor of &(d); and since 
the degree of S(d) is less than a &(d), (d) cannot be a factor of S(d). Thus, S(d) 
and &(d) are relatively prime. 

t This should be intuitively clear, since once 7 — (m — 1) bits are accepted in the 
ee the situation is one where the (remaining) input sequence is of a length 
ess than m. 
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S(d) such that 
U(d) + S(d) = T(d)#(d), 


where J(d) is some polynomial. In this situation the output has 
finite duration given by 7'(d), i.e., the periodic component of the solu- 
tion has been annihilated. To show this, we cite the following theorem.® 
Let U(d) and ® (d) be polynomials in GF(p). Then there are unique 
polynomials 7'(d) and S(d) inGF (p) such that U(d) + S(d) = T(d)®(d), 
where S(d) = 0 or S(d) is of lower degree then U(d). Recall that by 
suitably choosing initial conditions, S(d) can be any polynomial of 
degree m — 1 over GF (P). 

We turn to the important case of periodic inputs. The input se- 
quence U(d) can always be written in the form U(d) = P(d)/Q(d), 
where P(d) and Q(d) are relatively prime. Let the exponent of Q(d) 
be f, i.e., the period of the input is ¢. The d-transform of the output 
becomes 


Y(d) = [S@Q(@) + P(d)]/®@QM). (39) 


We consider first the case where ®(d) and Q(d) are relatively prime. 
If the numerator and denominator of eq. (89) are relatively prime, 
then, using the background material presented in Section VII, it is clear 
that the output is periodic with period N, where N = lem (1, p™ — 1). 
However, we will show that given P(d) and Q(d), there is a set of 
initial conditions for which 


S(@Q(d) + Pd) = Td)e(d), (40a) 


where 7'(d) has degree 1 — 1. When (40) holds, the output period is 1. 
Thus, assuming that all initial states are equiprobable, with probability 
p~™ the initial state will be such that the scrambler ‘‘locks up” and the 
output period equals the input period. (As we have previously men- 
tioned, this is a very undesirable situation.) To support (40) we cite 


the following theorem.’ There exist (unique) polynomials 7’(d) and 
S’(d) such that 


S'(d)Q@) + T'(d)e(d) = 1 (40b) 


only if Q(d) and #(d) are nonzero relatively prime polynomials over 
GF (p). Now multiply both sides of the above equation by — P(d) and 
let S(d) = — P(d)S‘(d) and T(d) = P(d)T’(d). Again we make use of 
the fact that S(d) spans the space of polynomials of degree m — 1 to 
guarantee that for every S’(d) there corresponds a S(d). We now 
summarize the above. 

Let the scrambler be defined by the primative polynomial 
1 — 3°”, c:d‘, and also suppose that the transform of the input to the 
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scrambler is P(d)/Q(d), where P(d) and Q(d) are relatively prime. It 
is also assumed that @(d) and Q(d) are relatively prime. For a particular 
set of initial conditions, the output period of the scrambler is the input 
period, J, where / is the exponent of Q(d). For all other initial conditions 
the output period is the least common multiple of 1 and p™ — 1. 

Our description is the same as Savage’s Theorem 1 with two dif- 
ferences—one superficial, the other crucial. Savage requires the poly- 
nomial h(d) = d™— >7,cd™* to be primative. However, 
1 — >°%,1 ed‘ and A(d) are reciprocal polynomials and, as we have seen, 
the reciprocals of primative polynomials are themselves primative with 
the same exponent. The second requirement is that ®(d) and Q(d) 
be relatively prime. This requirement, which is not part of Savage’s 
theorem, is essential for a complete description of scrambler behavior.* 

We will now show that the requirement that @(d) and Q(d) be rela- 
tively prime is satisfied whenever the exponent of Q(d) is not a multiple 
of p™ — 1. The proof is by contradiction. Suppose ®(d) and Q(d) are 
not relatively prime, then it is possible to write? 


Q(d) = R@)P(d) j= l,2,---, (41) 


where #(d) is a polynomial, with exponent 7, which is relatively prime 
to 6(d). The exponent of Q(d) is lem [7, p*(p”~) ], where k is such that 
pt <j S p*. Clearly, the exponent of Q(d) is a multiple of p” — 1, 
thus proving the desired result. Thus, when the input period is less 
than p”—" (the practical case), then @(d) and Q(d) are relatively prime. 
It is interesting to note that even if the input to the scrambler has an 
exponent which is a multiple of p” — 1, it may be that Q(d) and ®(d) 
are still relatively prime. For example, Q(d) can be the reciprocal 
polynomial to ®(d). 

Consider now the situation when Q(d) and ®(d) are not relatively 
prime. As above, we can then factor Q(d) in the form Q(d) =®/(d)R(a), 
j 2 1, where R(d) is either 1 or a polynomial relatively prime to 
@(d). From (39) and (41) the d-transform of the output is 


Y(d) = [S(d)#*(d)R(d) + P(d)]/P (AR). (42) 


Since by assumption P(d) is relatively prime with Q(d) = $4(d)R(a), 
the numerator and denominator of (42) are relatively prime. Since 
®*1(d) and R(d) are relatively prime, the output period is then the 
least common multiple of p*(p™ — 1) and r, with k given by p* <j 


“In other words, Savage states that, apart from the special case when the output 
period equals the input period, the output period is the lem (1, p”—!). This is not strictly 
true since, as we shall show, if @(d) and Q(d) are not relatively prime, the output 
period is not necessarily the lem(l, p—). 

t Since (d) is irreducible, we could not write Q(d) as a factor of &(d). 
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N 
Fig. 12—Cascade of N M-bit scramblers. 


+1<p*. Note that this result holds independently of initial 
conditions. * 

The above discussion provides a refinement of Savage’s basic result! 
by indicating that the output period is contingent on whether or not 
@(d) and Q(d) are relatively prime. It was shown that if the exponent 
of Q(d) is not a multiple of p” — 1, then &(d) and Q(d) are relatively 
prime. However, when Q(d) and #(d) are not relatively prime the 
output period must be determined from (42). 


IX. CASCADED SCRAMBLERS 


The cascade of identical scramblers provides an interesting example 
of when Q(d) and &(d) are not relatively prime. Suppose, as in I’ig. 12, 
we have n identical m-stage scramblers in tandem so that the output 
of the first is the input to the second and so on. Thus, assuming no 
lockup, the input to the second stage will have the same period as the 
free-running period of the second stage. With S(d) = 0 for all the 
scramblers, the output transform of the nth scrambler is 


Y(d) = U(d)/®*(d), (43) 


where U(d) is the transform of the input. Consider an example where 
U(d) = 1/(1+ d) = 1/Q(d) with p = 2. Note that the exponent of 
Q(d) is unity. The transform of the second output is 1/Q(d)#?(d), 
and we apply the results of the previous section to show that the output 
period is 2(2” — 1), 1.e., k = 1. Applying the above result to succes- 
sive stages produces the data in Table I, which shows the period of 
the output as a function of n for a binary scrambler (p = 2). Table I 
points out that adding a stage in cascade does not always double the 
output period. 

By considering each scrambler successively, we can comment on the 
output period of the cascade scrambler for arbitrary initial conditions 
and input, assuming Q(d) and ®(d) are relatively prime. Let the input 
to the first scrambler have exponent | < p”™ — 1. From Section VIII 


_ “If cancellation between the numerator and denominator in (42) were to occur, 
. (40) would imply that P = &(1 + SR®*—). Now since (41) states that Q = FR’, it 
is clear that P and Q have the common factor ®. This contradicts the assumption 
that P and Q are relatively prime. Thus, under the above conditions (i.e., Q and ® 
are not relatively prime), the initial condition cannot force the output period to equal 
the input period. 

t Results similar to ours were stated without proof in Ref. 9. 
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Table | 


n Output Period 
1 (2™ — 1) 
2 2 (2™ — 1) 
3 4 (2™ — 1) 
4 4 (2™ — 1) 
5 8 (2" — 1) 
6 8 (2 — 1) 
7 8 (2 — 1) 
8 8 (2™ — 1) 
9 16 (2™ — 1) 


we know that the probability of the output having period | is p-™. 
Otherwise, the output has period lem(l, p™ — 1). If the input to the 
second scrambler has period | we have the same situation as the first 
scrambler. However, if the output of the first scrambler has period 
lem (J, p™ — 1), the input polynomial to the second scrambler has de- 
nominator Q(d)®(d). By an argument analogous to that surrounding 
(42), it is clear that (40) cannot be satisfied, since Q(d)®(d) and ®(d) 
are not relatively prime; thus, the scrambler cannot. “lock up,” and 
applying the results in Table I indicates that the output of the second 
scrambler has period lem[1, p(p” — 1) ]. Thus, if a particular scrambler 
does not lock up, then no succeeding scrambler can lock up. The situa- 
tion is summarized in Table II for four binary scramblers in tandem. 
We assume in Table II that all initial states are equiprobable. 

We compare Table II to the serial scrambler in which all delay 
elements are combined into a scrambler that has 4m elements. 
With input period J, the output period is / with probability 2-4 and 
lem(l, 24" — 1) with probability 1 — 2-4”. Both the cascade and serial 
scramblers lock up and have period / with the same probability 
(2-4) ; however, since 


(t) the longest period of the cascade scrambler, Iem[/, 4(2" — 1)], 
is less than the largest period of the serial scrambler, 
[lem (J, 24" — 1)], and 


Table Il 
Output Period Probability 
l Q-4m 
lem(l, 2” — 1) Q-3m(]_ — Q-m) 
lem[l, 2(2" — 1)] 2Q-2m(] — 2-™) 
lem[l, 4(2" — 1)] 1 — 2-2" 
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(17) the probability of the cascade scrambler attaining its largest 
period (1 — 2-’”) is less than the probability of the serial 
scrambler attaining its largest period (1 — 2-4), 


the superiority of the serial over the cascaded scrambler in terms of 
spectral density is clear. 


X. PARALLEL SCRAMBLERS 


Serial scramblers have the property that if a single bit error is made 
in demodulation, then M errors will appear in the unscrambled output 
sequence, where M is the number of nonzero coefficients in the scram- 
bler primative polynomial. A parallel scrambler configuration has been 
proposed to ameliorate this error multiplication. In this section, we shall 
illustrate a spectral cancellation effect that can take place with parallel 
scrambling. For simplicity, we shall consider two parallel data streams. 
Suppose that the binary data, as in Fig. 13, are split into two data 
streams a, and b,, where an is the scrambled* version of a,, while 
bn = Gn_m@® b,. The an and the b, streams are then encoded for 
transmission over the channel. At the receiver, inverse operations re- 
cover the a, and the b, streams. A channel error in the a, stream will 
cause M errors in the a, stream and one errort in the b, stream. A 
channel error in the b, stream will cause a single error in the b, stream. 
Now suppose that a, and b, are Gray encoded so that a, is the least 
significant bit. The result is that the probability of error in the a, 
stream is much less than the probability of error in the b, stream; 
thus, the average number of errors in the a, and the b, streams will be 
decreased compared with serial scrambling. . 

Now we wish to examine the effect of ‘‘slaving’’ the b, stream to the 
a, stream. For our purposes it will be sufficient to code the scrambled 
output sequences into 1 and —1, ie., the transmitted data sequence 
is given by* (recall the notation of Section III) 


Cn = (242. — 1) + 7(2b, — 1) 
= 2an — 1 + j[2(An_m @ by) — 1). (44) 


As we have shown in Sections III and IV, both the line and envelope 
spectra critically depend on the discrete Fourier transform (pFT) of 
the cn, sequence, C(kQ). Unfortunately, it is not possible to express 


* At the inverse scrambler, a; is recovered as in the standard configuration, while 
bs is estimated as bn + Qn_m. 

T Since the estimated b. is formed as the “mod 2” sum of bn and Qn-m, & single 
channel error will only affect the 6, output once; however, the a, output will see 
the propagation of this error through the shift register. 

t The function (2a, — 1) maps “0”’ into ‘‘—1” and “1” into “1’’ and thus serves 
as a mapping from the scrambler output sequence to the transmitted (line) sequence. 
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Fig. 183—Parallel scrambler configuration. 


C'(kQ) directly in terms of the scrambler primative polynomial, which 
is a most difficult problem since, in the space of real-valued sequences, 
the prt is a linear operation and the scrambler output is a nonlinear 
function of the input [of course, the scrambler input and output are 
linearly related in GF (p) ]. However, in the present situation we are 
able to proceed since we are only interested in illustrating the possi- 
bility of spectral cancelling due to the parallel structure. We first 
indicate two simple relations between mod 2 operations and the cor- 
responding real variable operations: with a, b © GF (2), 


a@®b = (a — b)? (45a) 
a= a. (45b) 
Using (45) we write (44) as 
Cn = (24n — 1) + 9[2(Gn—-m — 2bD,dn—m + bn) — 1] 
= 2(dn + jdn—m) + 2jbn — Ajbpdn-m — (1+ 7). (46a) 


Let us consider the effect of the first term on the spectrum of c,. Now 
with LZ and C(kQ) denoting respectively the period of a, (and b,) and 
the prt of dn + j@n—m, we have* 


C(kQ) = [1 + je#0n!2) mk TA (KO) 
= 1 + efl(4/2)-@r/L) mk 4 (Q). (46b) 


Suppose that the scrambler produces a flat output spectrum (i.e., it 
would be a satisfactory scrambler if used solely in the serial mode), 
it is clear that C(kQ) will have periodically spaced nulls. Since the 
energy in the timing tone is given by 


Rss = |C(eQ) |2F (2) F (ws ~ bQ), (14) 


_ “For the purposes of our discussion, we, in effect, assume that b, = 0, 
1.€., Cn = On + 9Qn—m- 
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any amplitude tapering provided by |C(kQ) |? could impair the timing 
recovery system. From (46) we have 


|C(kQ) |? = [1 + sin(2r/L)mk} + [cos (2r/L) mk}? 
= 2[1 + sin (27/L)mk], (47) 


where k = 1, 2,---, LZ, andm = 1, 2, ---, M with M being the number 
of stages in the scrambler. The strength of the tone, described by (14), 
will be particularly attenuated when |C(kQ)| has a null at k = L/2, 
which corresponds to a frequency of w,/2. It is easy to see that for 
some values of ‘‘m,” the attenuation of the tone can be quite severe 
near w,/2 (i.e., k = L/2). 

In practice, the remedy is to change the value of m so that the null 
occurs as far away from w,/2 as is possible. Of course, this cannot be 
done prior to transmission since the exact value of w, is unknown. 

In this section, we have described a possible pitfall associated with 
the use of a parallel scrambler configuration. In practice, whether or 
not there is severe attenuation of the timing tone depends on the 
details of the pulse shaping and the operation of the phase-locked loop. 


XI. CONCLUSIONS 


We have examined several problems occurring in data-transmission 
systems that employ envelope-derived timing recovery, adaptive 
equalization, and self-synchronizing scramblers. Several conclusions 
have been reached regarding both the individual and joint action of 
these subsystems. 


(t) The performance of the envelope-derived timing recovery 
system can be significantly improved by narrow-zonal pre- 
filtering of the received signal prior to extracting the envelope. 

(it) The technique of ‘‘precessing” the data symbols in a phase-. 
modulated modem is sufficient to provide a timing tone in a 
large excess-bandwidth system, but does not provide a tone in 
a small excess-bandwidth system. 

(zit) A complete description was given of the output period of a 
cascaded scrambler as a function of the number of stages. Of 
interest are the facts that the output period does not neces- 
sarily double with the addition of a stage, and that if a partic- 
ular scrambler stage does not lock up, then no succeeding stage 
can lock up. 

(iv) It was demonstrated that the parallel scrambler configuration 
can, via spectral cancellation, cause the strength of the timing 
tone to be attenuated. 
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Analysis and Optimal Design of a Multiserver, 
Multiqueue System With Finite Waiting 
Space in Each Queue 


By B. A. WHITAKER 
(Manuscript received July 29, 1974) 


An analysis of a particular type of multiserver, multiqueue system 1s 
presented in which each queue has a finite number of waiting positions 
and the waiting positions are not vacated until service ts completed. Thus, 
several customers in one queue can be served simultaneously. The steady- 
state distribution of states is derived and ts used to obtain the probability 
of loss for each queue and the average delay of the system. This analysts 1s 
then used in the development of a design procedure to determine the 
minimum-cost configuration of waiting positions and servers to meet 
specified single-hour grade-of-service constraints. The results are applicable 
to the design of systems that utilize automatic call distributors. While this 
model does not include such effects as day-to-day variation and noncoinci- 
dence of peak loads among trunk groups, nevertheless the results properly 
reflect for the first tume the interactions among the trunk groups terminated 
on the automatic call distributor and the attendants at the automatic call 
distributor. 


1. INTRODUCTION 


The purpose of this paper is to present an analysis of a particular 
type of multiserver, multiqueue system in which each queue has a 
finite number of waiting positions and the waiting positions are not 
vacated until service is completed. In particular, the steady-state dis- 
tribution of states are derived, and expressions for the probability of 
loss for each queue and the average delay are given. It is shown that 
the queuing model described is representative of systems characterized 
by a finite number of trunk groups that carry calls to a group of 
attendants who then perform some service for the caller. (One such 
system is used in the directory assistance service provided by the 
telephone company.) Results are given to illustrate the effects of vary- 
ing the number of servers and number of positions in each queue. 
Finally, a design procedure to determine the minimum-cost configura- 
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tions for such systems under various grade-of-service constraints is 
developed. This procedure ignores such effects as day-to-day varia- 
tion, noncoincidence of peak loads among incoming trunk groups, and 
retrials of blocked calls, which should be investigated in the develop- 
ment of procedures for traffic engineering and administration. However, 
as in most cases, it is difficult, if not impossible, to obtain analytical 
results with these effects included. This paper should provide useful 
insight that can later be incorporated in a complete traffic engineering 
procedure. 

The system analyzed consists of J input queues each with a finite 
number of waiting positions, NV; (¢ = 1, ---, l), which have full access 
to M servers. When an arrival seizes one of the M servers, it does not 
vacate a waiting position, but remains in the position until its service 
has been completed. This characteristic, which allows calls that are 
not at the head of a queue to be in service, distinguishes this system 
from the usual queuing system. In particular, the system is no longer 
completely described by the number of calls in each queue since a 
record of the number of calls from each queue that are in service must 
be kept. An arrival that finds all the waiting positions for its queue 
occupied is cleared or lost from the system. An arrival finding no idle 
servers but at least one vacant waiting position in its queue enters the 
queue and is delayed until its service begins. In the context of directory 
assistance systems, the input queues are the trunk groups and the 
waiting positions are represented by the trunks. The information 
operators are the servers. 

In telephone traffic theory, the described system has been referred 
to as a combined loss-and-delay system. Previous work in this subject 
can be segmented into three parts: 


(t) One input flow—one queue.!-® 
(tz) Several input flows—one queue.°—"! 
(zi7) Several input flows—several queues.!2—!4 


The last segment, of which this analysis is a part, has been investigated 
by Kithn. He analyzed systems with g > 1 queues, each with a finite 
number of waiting positions s; (¢ = 1, 2, ---, g). Associated with each 
queue is a Poisson arrival process with mean rate \;, which is assumed 
to be independent of the others. An arrival that finds all waiting 
positions in its queue occupied is lost. Arrivals that are not cleared from 
the system are served by one of n servers. The service time distribution 
for the ith server is exponential with a mean rate e;. When a server 
becomes idle, queue 7 is chosen to receive service with probability p;. 
Within a queue, calls can be selected randomly, first-come, first-served, 
or according to a priority scheme. 
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As Kihn indicates, analytical solutions to these systems exist only 
for special cases. In general, the linear equations representing the 
equilibrium conditions must be solved numerically for the particular 
values of the parameters. However, Kiihn gives a solution to one 
particular system, which will now be discussed. In this system, it is 
assumed that the service rate for all servers is identical (e; = «) and 
the interqueue discipline is defined by the p,’s that are 


Z; 

g 
D 2 
k=1 





Pi = (PSd, 989); 


where Z; is the number of waiting positions occupied in the jth queue. 

The system analyzed in this paper is an extension of the one ex- 
amined by Kihn, since an arrival does not release a waiting position 
until his service has been completed. This complicates the state 
analysis, since it is now not sufficient to know the number of servers 
that are busy to determine the equilibrium equations; information as 
to the number of calls from each queue that are in service must be 
included. 

Kiihn indicates also that, for the above interqueue discipline, the 
waiting time distribution can be found numerically only for small 
systems. The calculation of the waiting time is complicated by the 
fact that an arrival’s waiting time is influenced by the number of 
arrivals that occur after it has entered the system. More is said about 
this difficulty later. 


Il. MATHEMATICAL FORMULATION 


In this section, a mathematical model of the queuing system is 
formulated. Equilibrium equations are given and their solutions 
derived. 


2.1 Queuing model 


The queuing system consists of J input queues, each with a finite 
length denoted by N;, 7 = 1, 2, ---, 1. Requests for service arrive at 
queue 7 according to a Poisson distribution with mean rate, \;. If we 
let A,(¢) denote the number of arrivals at queue 7 in (0, é), then 


(At)* 4, 
P[A;(#) = k] => as dit (k = 0, 1, 2, oe -), (1) 
The arrival process at queue 7 is assumed to be independent of the 
arrival process at each of the other queues. Arrivals from each queue 


have full access to a group of M servers. The service time distributions 
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of the servers, denoted by A(t), are independent and identical ex- 
ponential distributions with mean service rate y, i.e., 


1 — e#! (¢ = 0) 
HOG (t <0). @) 
Since the arrival process is Poisson and the service process is ex- 
ponential, the queuing model is a multidimensional birth-death process. 

Arrivals to queue 7 that find NV; waiting positions in queue 7 occupied 
are not allowed to enter the system. If an arrival to queue 7 finds at 
least one unoccupied position in queue 7 but all M servers busy, it 
enters the queue and waits as long as necessary for service. Within 
queue 7, arrivals enter service on a first-in, first-out (FIFO) basis. An 
arrival that finds at least one unoccupied position in its queue and at 
least one free server immediately enters service. When an arrival 
enters service, it does not release a waiting position but remains in 
the queue until its service has been completed. Hence, the word 
“queue” is being used in a nonstandard manner and refers to the 
number of calls waiting for service and in service. As discussed earlier, 
an example of such a queue is a group of trunks that carry calls into 
a switchboard. 

The interqueue service discipline—the order in which the queues 
receive service—is characterized by the number of calls waiting for 
service. When a server becomes free, queue 7 receives service with a 
probability, p:, which is the proportion of queue-7 calls waiting for 
service. If we denote the number of calls in queue 7 by n; and the 
number of calls in queue 7 that are in service by m:, then this prob- 
ability, which is dependent on (m1, m2, +--+, M1, M1, °**, mi) = (n, m), 
can be expressed as 





2 eas : 2 fo . U 
nN; Mi _ Ni Mm; (x nj = M +1) 


(i = 1, 2, a8 acl). (3) 


The effects of other interqueue service disciplines have been in- 
vestigated, but will not be discussed here. 

The fact that arrivals remain in the queue during service dis- 
tinguishes this queuing system from the standard system, since the 
amount of information required to fully describe a state of the system 
is increased. The system is also complicated by the fact that the inter- 
queue service discipline is state-dependent. However, as is shown in 
later sections, this “complication” leads to a closed-form solution of the 
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equilibrium equations, which is generally not the case for such systems. 
The equations of equilibrium are given in the following section. 


2.2 Equilibrium equations 


The system described in the previous section is characterized by a 
finite number of states that indicate the number of calls in each queue 
and, of these calls, the number receiving service. We denote by (n1, ne, 

++, M1, M1, +++, M1) the state in which there are n; calls in queue 7 
with m; of these calls in service. For notational simplification, we also 
refer to the state in vector notation as (n, m). In this notation, (n:;, m) 
represents the state (mi, m2, --:, mi +1, ++, mi, mi, *-+*, mi) and 
similarly (n;_,m) = (n:, ---, ni — 1, ---, mi). It should be clear that, 
if the total number of calls in the system is less than or equal to M, 
then n; = m; for all 7. 

Assuming stationarity, let P(n, m) be the probability that at an 
arbitrary instant of time the system is in state (n, m). Moreover, if the 
arrival processes to the system are Poisson, the equilibrium-state 
distribution {P(n, m)} at an arbitrary instant is equal to the equilib-. 
rium-state distribution at the instant of an arrival. By equating the 
rate into a state to the rate out of a state, we can write the equilibrium- 
state equations where we have introduced the function 


1 x>0d0 
UD=j19 <0 


to include the boundary conditions and a; = \;/u: 
l 
| >» Law(Ni — ni) + nal P(n, m) 
i=1 
U U 
=; > aiu(n:)P(ni-, mi) + Xu (ni + Ll)u(N;i — ni)P Maz, mix) 


(<a) (0S n; S Ni) @4=1,2,:-:,t (4) 
| y Law(N; — ni) + mi} P(a, m) 
= = awu(n;)P (ni, m,;_) + E miu = ni) P(ni+, m) 


+ ¥ (m+ Du(Wi — ndu(m)P (is, mis) 


° 
HI 
ran 

oe 


(x =m) OSweNS. €H Se) 
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Ps * Law(Ni — ni) + mal| P(n, m) = > awu(ni)P(ni_, m) 


cn u(N; — n:)P(niz, m) 


Eo Gat 1)(m; + 1 — m,) 
pp» (= = 
pHi Y~mtl—- M 
ya> 


u(N; — niu(m,)P(nizy, Miz 5) 


1) (O<n5N) i=1,2,---,h (6) 


Equations (4) to (6), together with the normalization condition 


Ni Ni min (M,N1) min (M,N1) 
pap 2 oD vee DD Pia, ia) = 4, (7) 
m =0 m=0 m,=0 mi=0 
where all nonexistent states, such as states in which both n; = 0 and 
m; > 0, in the sum are assumed to have probability zero, determine the 
equilibrium-state distribution. 


2.3 Steady-state solution 


Since the process described by the equilibrium equations is a 
finite-state birth-death process in which the arrival rate into the 
system is always less than the service rate of the system as a result of 
overflow from the finite queue, a unique solution to eqs. (4) to (7) 
exists, and the solution is a genuine probability distribution. This 
solution is given in terms of P(0, 0) by 





l ni U 
I P@,0) (2mm) (8) 
t=1 Ni: i=1 
aa U l 
P(n, m) = Ya M I av 
i=1 1=1 
Ni — M1, + °°, M1 —~ My (M1, Ma, ++ *, M1 


M! M21-™ P(0, 0) 


(dm > Mt), (9) 


where P(0, 0) is determined from (7). The general solution was de- 
termined from examination of various small systems. By substitution, 
it can be shown that this solution in fact satisfies the equilibrium- 
state equations (4) to (6). 

When the number of calls in the system is less than or equal to the 
number of servers, no one is waiting for service and the queues have no 
interaction. This fact is shown in (8) by the product form of the solu- 
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tion. However, as the number of calls increases to the point where calls 
are waiting in more than one queue, the queues no longer are acting 
independently. 

Since the queues behave independently as long as there are free 
servers, we would expect that, as the number of servers was increased, 
the system would approach / independent loss systems. By examining 
the marginal probabilities, it can be shown that this is, in fact, true 
when M = >°}_, N;. That is, the marginal state probability of mi 
calls in queue 1 is 


Pom) = E+E Pla, m) = yi, (10) 
> ak/k! 





which is identical to the state probability for a pure loss system. 

If no calls were blocked from the system, then the system would act 
as a pure delay system with the offered load a = >-}_, a;. This is 
easily shown by taking the limit of (8) and (9) as N;—> © for all 
$=1,2-+1. 

We first consider the case in which }7/_, n; S M. Since the number 
of calls in each queue is unrestricted, it is easily seen that the multi- 
nominal expansion of (a1 + --- + a1)" divided by (S~j-; 7.) ! can be 
obtained from (8). That is, 


1 
More, Bete Pla, m) = FOIE peo) (yu <M): 
Hence, 


l ak 
P (> fie k) = 2 P00) (<M). (11) 


For the case in which >(}_, n; 2 M, first note that, by the Vander- 
monde convolution of multinomial coefficients, 


Ut 
> ‘ > Ni ar M M 
or » i=1 
mi mt 
N1— M1, *°°, Ni — M1 mi, °° *, M1 
rf 
De ni 
= er - (12) 
n1, N2, ee NI 
Therefore, 


(2%)! i ai" TPO, 0) 


P(n) = 2 P(a, m) = Mt agene 


(Reem), as 
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and thus we can denote this state probability by 


1 L af 
P( Yaak) =, » Saran Tn , P(0, 0) 


2 n= 
(k 2M). (14) 
Consequently, 


l 
we P (x n= k) = THER a P(0, 0) (k 2M). (15) 
i=1,- 


The normalization constant is then expressed as 


M a*® —1 
Po =| Da i+ Sy arare | 0 <a<M. (16) 


Hence, comparison of (11), (15), and (16) with the pure delay system 
completes the proof. 


2.4 Blocking and delay probabilities 


In the analysis and design of queuing systems, performance mea- 
sures for each configuration must be calculated. In telephone traffic 
theory, these performance measures are generally referred to as 
“grades of service.’”? Two such measures of the grade of service are: 


(4) For loss systems, the blocking probability or probability of loss. 
(ii) For delay systems, the average delay experienced by calls that 
enter the system. The average delay W(s, a) for a pure delay 
system is expressed in terms of the Erlang delay formula as 


C(s, a) | 
(s — @)p 


In the system described in Section 2.1, the blocking probabilities 
for each input queue, the average delay experienced by calls that enter 
the system, and the average delay of only those customers who experi- 
ence a, positive delay are important characteristics to be examined. The 
latter is not used in the remaining analysis. 

The blocking probability for queue 7 is defined as the probability 
that an arrival to queue 7 finds NV; calls in the queue. This probability, 
which is denoted by B:(N, M, a), is a function not only of the number 
of positions in queue 7 and the offered load to the queue, but also of 
the traffic load offered to each of the other queues, the number of 
positions in each of these queues, and the number of servers. Recalling 
that, for systems with Poisson input, the state probabilities at an 
arbitrary instant are equal to the state probabilities at arrival times, 


W(s, a) = (17) 
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B;(N, M, a) is calculated from the marginal distribution for queue 
1 as 


BAN, M, a) 
= ey wee Dee DB Pm, +++, Ni, +++, m1, m). (18) 
ni Ni-k nit1 mi ml 


The average delay experienced by calls that enter the system (suc- 
cessfully occupy a waiting position in their queue) is denoted by 
D(N, M, a). The delay calculated for this system is the overall mean 
waiting time measured from an arrival’s entry into its queue until its 
service begins. Hence, it does not include service time of the call. It 
should also be noted that arrivals into the system that find at least one 
free server experience no delay. 

Since the average number of calls waiting for service must be finite 
and the mean waiting time is finite as a result of the loss structure of 
the system, the well-known equation of Little,° L = \W, can be used 
to calculate D(N, M, a). In particular, we must define our ‘“‘queue 
length” as the total number of calls waiting for service and ‘“‘d’’ is 
defined as the effective arrival rate into the system. Hence, 


a > see x] (5 a M1) PCa, m) 

ZN; 

DN, M, a) = ———_, -—_: (19) 
u AL] — BN, M, a) ] 


Calls that are blocked from the system do not enter the queue and 
hence do not affect the average queue length. Consequently, they are 
not included in the arrival rate into the system. The numerator of (19) 
is the average number of waiting calls. It should be apparent that the 
average delay for any particular queue can easily be obtained by using 
the appropriate marginal probabilities. Also, the conditional average 
delay, the average delay experienced by only those that must wait, 
is found by dividing (19) by the probability of being delayed. 

For some design purposes, it might be deemed necessary to constrain 
the probability of waiting longer than some time, to, to be less than a 
specified value. In this case, the waiting time distribution for each 
queue must be obtained. For the interqueue discipline examined for 
this system, the calculation of the waiting-time distribution is ex- 
tremely difficult. (Kiihn™ mentions that, for his problem, numerical 
techniques can be used for very small systems, after which approxi- 
mate methods must be formulated.) The difficulty in determining the 
waiting-time distribution lies in the fact that the time a particular call 
must wait for service is not just a function of the number of calls in 
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the system when it arrives, as is usually the case. In particular, the 
waiting time of a call is related to the number of calls that arrive after 
the particular call enters the system, since queues are chosen for 
service according to their queue lengths at a service completion, so 
that later calls can ‘“‘pass” earlier ones. For this reason, the waiting 
time distributions are not calculated in this study. 


2.5 Macrostate analysis 


As discussed in the previous section, blocking probabilities and 
average delay values are of interest to system designers. Using the 
state probabilities given in (8) and (9), it is possible to obtain not only 
the overall average delay, D, as shown in (19), but also the average 
delay, D;, for queue 7. In certain types of design, we might want to 
engineer the system so that the average delay in every queue is less 
than a specified level. In such cases, D; would be needed. However, for 
this study, we consider only the overall average delay. 

Therefore, it is apparent from (19) that, for computational purposes, 
we only need to know the number of calls in each queue without dis- 
tinguishing between those in service and those waiting. If we denote 
by n the state (m1, m2, ---, m1), we can find the steady-state prob- 
abilities P(n) from (8) and (9). Of course, we could have written the 
state equations directly and solved this easier set of equations.!” 
However, for further studies, it is essential to know the probabilities 
P(n, m). 

Since the state probabilities for (n, m) in which }°}_,; S M are 
independent of m, P(n) = P(n, m). To obtain P(n) for Di-in: > M, 
we sum P(n, m) given by (9) over all possible values of m. Using the 
Vandermonde multinomial convolution, we find that 


(3m)! 


Pea) = SE Pam) = aes ae I PO), 20) 


where P(0) is the normalization constant. We can calculate the block- 
ing probabilities and the average delay as before. The number of 
states is [[$-, (Vi + 1). 

It should be repeated that the macrostate probabilities are of use 
only if one is interested in the overall mean delay. To calculate the 
individual average delays, one must use the microstate probabilities. 


Ill. RESULTS 


In this section, we investigate the effects of varying N; and M on 
the blocking probabilities for each queue and the overall average delay 


604 THE BELL SYSTEM TECHNICAL JOURNAL, MARCH 1975 


of a call. In particular, these effects are illustrated by comparing the 
results obtained from the analysis presented in this paper, for a partic- 
ular example, with the results obtained if independence is assumed. 
This ‘independence assumption” is often used in practice to deter- 
mine the number of positions for each queue and the number of 
servers needed. In essence, this assumption permits a designer to 
design each queue (trunk group) independently of the other queues 
and the number of servers, and the number of servers is determined 
assuming that no calls are blocked in the queues. It is shown that this 
assumption is generally not even a good approximation. Finally, four 
properties that are used in a design procedure established in the next 
section are postulated. 


3.1 Comparison of results with independence assumption 


In the engineering of automatic call distributor systems, a traffic 
engineer generally dimensions each trunk group using the Erlang loss 
formula (assuming that it is independent of the other groups and the 
number of attendants) and often determines the number of attendants 
required from the Erlang delay formulas using the total offered load 
to the queues (assuming no blocking in the queues). This procedure is 
clearly invalid, but up to now an exact procedure has not been avail- 
able. As a means of illustrating the significance of the results presented 
in this paper, we now compare, for a particular system configuration, 
the system characteristics that a traffic engineer would expect to 
obtain using the independence assumption and what he really will 
find. Of course, the interqueue service discipline will affect these 
results in a way that will be described in later work. The actual opera- 
tion of such systems is quite complicated, and is not readily character- 
ized by any of the disciplines usually used such as first-in, first-out, 
random, and last-in, first-out. However, the results of simulations 
indicate that the discipline presented here is a good approximation of 
the actual method of operation. 

For purposes of this comparison, we assume a simple system with only 
two queues: the first with an offered busy-hour load of 10 erlangs; the 
second, 5 erlangs. It is further assumed that the queues have coincident 
busy hours and that the average holding time per call is 30 seconds. 
Assuming that a P.01 grade-of-service constraint has been placed on 
each group, the number of trunks required, if independence is assumed, 
would be N; = 18 and Nz = 11. (These numbers can be found from 
tables of the Erlang B formula.) If it is then required that, on the 
average, no call must wait longer than 3 seconds for an answer, we 
find, from the Erlang delay formula, that M7 = 19 (assuming no 
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blocking in the groups). The traffic engineer would expect this system 
to have the following characteristics: 


Blocking probability for group 1 = 0.0071. 
Blocking probability for group 2 = 0.0083. 
Overall average delay = 1.83 seconds. 


However, analyzing this system configuration with the results of this 
paper, we find that the service levels (which stated above are approxi- 
mations to the actual levels) would be 


Blocking probability for group 1 = 0.014. 
Blocking probability for group 2 = 0.018. 
Overall average delay = 0.949 second. 


The directions of the changes are intuitively obvious, since longer 
holding times in the queues result in more blocked calls and the higher 
blocking levels decrease the load to the servers, which in turn results 
in a lower average delay. It should be noted that the trunk groups 
are performing at unsatisfactory levels, but the overall average delay 
has been decreased and is considerably under the required level. Often, 
customers who have such systems measure only the delay or speed 
of answer and periodically remove attendants if they feel that the speed 
of answer is not above the required level. Unfortunately, such a 
customer generally does not realize the effect of removing attendants 
on the blocking probabilities on his incoming trunks and consequently 
on other network customers. 

As an example of customer behavior, consider the system discussed 
above. The customer, having measured the average delay and finding 
it to be considerably under his required level, would most likely re- 
move two attendants. The average delay would then become 3.01 
seconds, but the blocking probabilities increase to 0.028 for group 1 
and 0.023 for group 2. Hence, even though the delay constraint is 
essentially satisfied, the probabilities of loss are more than double their 
desired levels. 

If, instead of the P.01 service level, P.05 or P.10 had been chosen. 
for the above delay constraint, the independence assumption would 
generally give a configuration that would satisfy all service con- 
straints. The reason for this is that the higher blocking levels decrease 
the offered load to the attendants, and consequently a very small delay 
results. This delay is small enough that it has little effect on the hold- 
ing time of the calls and, hence, the offered load to queue 7 is ap- 
proximately a;. Therefore, the resulting blocking probability, although 
larger than the Erlang loss probability, is generally in the acceptable 
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range. However, the configuration would not be optimal because too 
many attendants are provided. In summary, as the blocking prob- 
abilities increase, the discrepancies between the Erlang loss formula and 
the formula given by (18) decrease, but the discrepancies between 
average delays increase. 

Several variations in the procedures to engineer these systems exist 
in practice. Generally, in determining the number of attendants, the 
measured offered load into the attendants is used. This load accounts 
for the blocking in each trunk group. For the above example, the 
measured offered load would be 14.79 erlangs (if we assume that the 
assumptions made in this analysis are valid) and, from the Erlang 
delay formula, an average delay of 1.58 seconds results, which is still 
higher than the actual delay. The reason for the discrepancies is that 
the offered load to the attendants is no longer Poisson as a result of 
the blocking in the groups. In fact, the variance of this offered load will 
be lower than that of the Poisson load, since the peakedness of the 
traffic has been decreased by clipping. Hence, since the actual average 
offered load and variance of the load are lower, this leads us to postulate 
the following property : 


Property 1: 


7 7 l 
W(M, a) 2 D(N, M,a) where a= > ai. 


t=1 


This property is illustrated in Fig. 1. The equality holds in the limit 
as N; > for all 7 = 1, ---, J, as shown previously. The significance 
of this property is that we now have a method of obtaining an upper 
bound on the average delay for the combined system. 

Another variation that is sometimes used is to add the speed of 
answer into the offered load to each group. If we add the average 
delay of 0.949 second to each call and use this new offered load in the 
Erlang loss formula, the blocking probabilities that result are 0.0091 
and 0.01, for groups 1 and 2, which are still lower than the actual 
blocking. 

The discrepancies result because the Erlang loss formula assumes 
exponential holding times on the trunks but, in fact, the holding time 
for the combined system is the sum of an exponential distribution and 
the delay distribution. Also, the holding times of calls in the system 
are no longer independent (unless >°j., ni < M). The variance of this 
new service time distribution is higher than that of the exponential 
service time distribution and, of course, the mean is larger. Therefore, 
one would expect the average queue length to be larger which, in 
turn, implies an increase in blocking. 
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Fig. 1—Average delay as a function of number of positions in one queue. 


With these facts in mind and as a result of empirical evidence, we 
postulate a second useful property : 


Property 2: 
B(Ni, ai) S BN, M, a) i= 1, hel), 


where B(N;, a:) is the Erlang loss formula for N; servers and an offered 
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Fig. 2—Blocking probability as a function of number of positions in the queue. 


load, ai. The equality holds when M = ZN, since the queues then 
behave independently (there is no delay, so in fact the holding time 
per call is exponential). This property is illustrated in Fig. 2. In- 
tuitively, one would expect B;(N, M, a) to be larger since, as a result 
of a positive delay added to each call, calls hold the trunks longer, 
therefore increasing the probability of an arriving call finding all 
trunks busy. The significance of Property 2 is that a lower bound on . 
the number of trunks required for a given service level can be found 
using the Erlang loss formula. 
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3.2 Effects of varying N, and M 


For a given input process and service time distribution, the designer 
can affect the blocking probabilities for a queue or the average delay 
by changing the number of positions, N;, in a queue and/or changing 
the number of servers, 1. We first investigate the effects of varying 
the number of servers, /. As noted earlier, when M has been increased 
to the point where >"}.,N; = M, the average delay becomes zero 
and the queues behave independently. The blocking probability for 
each queue is then given by B(N;, a:) which, by Property 2, is a lower 
bound on the blocking for any value of M. What is of importance, 
however, is: Do the blocking probabilities and the average delay 
monotonically decrease to their lower bounds as we increase M to the 
value >°}., N;? Empirical evidence, such as shown in Figs. 3 and 4, 
indicates that the answer to this question is yes. We postulate this 
property as: 


Property 3: 


BAN, M + 1, a) < B;(N, M, a) (@=1,---,) 
D(N, M +1, a) S D(N, M, a). 


Intuitively, one would not expect that increasing the number of 
servers in a system would increase the average delay. Moreover, a 
decrease in the holding time of calls would imply that the average queue 
lengths would decrease, which would result in a decrease in the blocking 
probability for that queue. However, this decrease in blocking results 
in an increase in offered load to the servers but, as we postulate, this 
increase is less than the marginal carrying capacity of the added server. 
The significance of this property is that, with added servers, not only 
is the average delay decreased but also the blocking probabilities are 
decreased ; that is, adding servers improves the service performance of 
the servers and of the queues (trunk groups). 

The other system parameters that may be varied to improve system 
performance are the numbers of positions in each queue. Supported by 
quantitative evidence, such as given in Figs. 1 and 5, and by intuition 
we postulate the following: 


Property 4: 
BAN.,, M, a) S BN, M, a) 
B(Niz, M, a) = BN, M, a) j# 1. 
DNs, M, a) = DIN, M, a) 


The first part of this property states that, if we increase the number of 
positions for calls to occupy in a given queue, then the probability of 
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Fig. 3—Average delay as a function of number of servers. 


loss for that queue is decreased. (This, of course, is true in pure loss 
systems.) The intuitive argument is that calls that previously found 
N; calls in queue 7 were blocked, but now are not. Therefore, the 
number of calls blocked is decreased. However, as the third part of 
this property implies, the average delay of calls is increased as a 
result of this increase in calls from queue 7. The intuitive counter- 
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Fig. 4—Blocking probability as a function of number of servers. 


argument is that this increase in holding time per call might result in 
an increase in blocking for queue 7. But we postulate that the increase 
in delay is not substantial enough to eliminate the increased efficiency 
obtained in queue 7 by the addition of a position. 

However, for the other queues, the number of positions remains 
fixed, and this increase in average delay results in a larger traffic level 
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Fig. 5—Blocking probability in one queue as a function of number of positions in 
another. 


per position being placed on the queues. Hence, this increase in load 
causes the average queue length to increase and, consequently, causes 
the blocking probabilities to increase. This fact is stated in the second 
part of Property 4 and is illustrated by the example given in Fig. 5. 
Equality in the three statements holds only in the limit as N; goes to 
© for allt = 1, 2, ---, 1. 

The significance of Property 4 is that the loss probability for a 
given queue cannot be reduced by adding positions to any other queue. 
Therefore, by Properties 3 and 4, we see that the loss probability for a 
given queue can be reduced only by increasing the number of servers 
or by increasing the number of positions in that queue. 

The four properties postulated indicate relationships between the 
system parameters and the system characteristics. Proofs for the 
simplest cases (i.e., ] = 1 for all the properties except for the second 
part of Property 4 for which | = 2) are given in Ref. 18. The proofs 
for the general cases have not been constructed because of the dif- 
ficulties involved (e.g., the proof of the second part of Property 4 
required 17 pages for 1 = 2). However, based on the intuitive explana- 
tions given, empirical evidence, and these proofs, I feel that the 
properties are valid in the general case of 1 queues. To obtain an optimal 
configuration (minimum cost), one must balance the cost of servers 
against the cost of positions for the queues in such a way that all 
required service levels are met. These properties are used in the next 
section in the development of a procedure to determine this optimal 
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configuration. The algorithm and proof of convergence given below 
assume the validity of these properties. 


IV. SYSTEM DESIGN WITH GRADE-OF-SERVICE CONSTRAINTS 
4.1 The design problem 


In determining the optimal design of a queuing system, one is 
generally interested in minimizing the operating costs in such a way 
that specified grade-of-service constraints are met. These constraints 
are a function of the particular queuing system under study. In the 
pure loss system, the grade of service is measured by the probability that 
a call is lost or blocked. Hence, the optimal system configuration is the 
minimum number of servers that satisfies the constraint, B(s, a) S 6. 

In delay systems, at least three measures of service are useful. The 
first is the average delay experienced by a call in obtaining service. 
The second measure is the ‘‘extremal’’ delay—the probability that the 
delay for any call exceeds a specified limit. Finally, the average delay 
of only those customers who experience some delay is a useful measure. 

However, in the combined loss and delay systems described in Sec- 
tion II, the determination of an optimal configuration is not as straight- 
forward. We will measure the grade of service of the system by 


(t) The blocking probability for each queue. 
(ii) The average delay of all calls that enter the system. 


The blocking probability for a particular queue is dependent on the 
number of positions in each queue and the number of servers, and the 
average delay depends on these same variables. A procedure must be 
developed to balance these measures of congestion in such a manner 
that the costs of the system are minimized. 

More formally, the problem can be expressed as the following 
nonlinear problem in integer programming. We denote the monthly 
cost of a waiting position in queue z by C; and the monthly cost of 
each server by C. It is assumed that C; and C are positive, finite 
numbers. The blocking objective for queue 7 will be denoted by 6; and 
the average delay objective by d. The following assumptions have been 
made: The system is engineered for the system busy-hour traffic load 
and the busy hours for the queues are coincident. The optimization 
problem is then expressed as: 


Minimize the cost 


Tt 
Z=>DCNi+ CM 
il 
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subject to 
(I) BN, M, a) = b; (2 cose Oe eg l) (21) 


D(N, M,a) <d (22) 
N.(t = 1, 2, ---, 0, M = 0 integers. 


Since no efficient algorithm to solve a general nonlinear integer- 
programming problem exists, a procedure was developed that utilizes 
‘ the four properties stated in Section III. 


4.2 Optimization procedure 


In this section, a procedure is developed that determines an optimal 
solution to the nonlinear integer-programming problem expressed by 
(I). The procedure is a direct-search routine based primarily on the 
properties presented in the previous section. A description of the 
procedure is now given and is followed by a concise summary of the 
algorithm. Figure 6 is a flowchart of the algorithm. 

The first step in the algorithm, as in most mathematical program- 
ming algorithms, is the determination of a feasible solution to the 
problem. To obtain an initial feasible solution to (1), we utilize Proper- 
ties 2 and 3 of the previous section. In particular, by Property 2, we 
know that the minimum number of waiting positions for queue 7 can 
be determined from the Erlang loss formula, which is easily computed 
from a recurrence relationship.!® We begin the search for an initial 
feasible solution with 


Ni = min {n|B(n, ai) S bi}, 


since it has been shown that, in fact, a feasible solution to (1) exists. 
That is, (N®, M2), where I = S}., N, is a feasible solution since 
D(N®, M,a) = Oand Bi(N®, M, a) = B(N, a) fori = 1, 2, ---, 1. 
However, since this solution will generally not be near the optimal 
solution, the search will not begin at MZ but instead with Mo), which 
is the minimum value of M that satisfies the constraint: 


W(M, a) $d. (23) 


Mo) is the number of servers that would be selected if the Erlang 
delay formula with an offered load a = )°}., a; were used. By Property 
1, it is seen that if (23) is satisfied, then (22) will also be satisfied. 

Using the set of parameters (N, Mo), the system characteristics, 
{B,(N, Mo), a)} and D(N™, Mo), a) are determined. One of two 
results occurs: 
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(t) The parameters (N®, M o)) satisfy the constraints of (I), in 
which case an initial feasible solution has been determined. We then 
proceed to find the feasible solution of minimum cost. 

(ti) At least one of the blocking constraints (21) is violated. (As 
noted above, the delay constraint will be satisfied.) By Property 3, 
we know that the addition of a server will reduce the blocking prob- 
ability for each queue and that, by the addition of enough servers, a 
feasible solution can be obtained. We denote this initial feasible solu- 
tion by (N, M) and note two interesting properties of this solution. 


(a) By Property 2, N{ is the minimum number of positions that 
must be considered for queue 7. 

(b) M® is the maximum number of servers that must be con- 
sidered. This is true since a further increase in the number of 
servers can be justified only if some waiting positions can be 
eliminated, and the positions are already at their minimum 
levels, N. 


The next step in the algorithm is to attempt to improve the initial 
feasible solution. Since by construction we are initially at the maxi- 
mum number of servers, we attempt to decrease the costs by decreas- 
ing the number of servers while maintaining feasiblity. To maintain 
feasibility, it may be necessary to add waiting positions to certain 
queues. If a feasible solution is found, then it will be an improvement 
only if the accumulated cost of those servers removed (since the last 
feasible solution) is greater than the accumulated cost of all waiting 
positions that have been added to maintain feasibility. Hence, as we 
remove servers, one of three things results. 


(z) All the constraints of (I) are satisfied. If the accumulated cost 
of removed servers is greater than the cost of all waiting positions that 
have been added, this new feasible solution represents a cost improve- 
ment and should be stored as the tentative “‘optimal” solution. The 
accumulated costs are set to zero, a server is removed, and the search 
continues. If the cost of servers is less than the cost of positions, then 
we reduce the number of servers by one and continue the search for a 
solution with lower cost. 

(it) The delay constraint (22) is violated. In this case, we stop the 
search and the tentative optimal solution is the global* optimum. (A 
justification for stopping the procedure is given later.) 

(iz) At least one of the blocking constraints is violated. In this case, 
we add one position to each queue in which the corresponding block- 


“I have taken the liberty of using “global’’ since, in fact, the procedure does 
produce the global optimum 7f the four properties are true in the general case. 
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Fig. 6—Design procedure. 


ing constraint is violated and increase the accumulated cost of added 
positions appropriately. If the cost of the additional positions is less 
than the accumulated cost of the servers that have been removed, we 
determine if this solution is feasible. If it is feasible, we proceed as in 
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(z). If it is not feasible, either (iz) or (¢77) must be true. If the cost of 
the positions is more than the cost of servers, then no feasible solution 
with lower cost can be obtained with this number of servers and, hence, 
we reduce the number of servers by one and calculate the system char- 
acteristics with the new set of system parameters. Then either (2), 
(ii), or (717) must be true, and the appropriate action is then taken. 


In summary, the procedure removes servers until either the delay 
constraint is violated, in which case it terminates, or a blocking con- 
straint is violated. If a blocking constraint is violated, waiting posi- 
tions are added, and an additional server may be removed, depending 
on the incremental costs and on the feasibility of the tentative solution. 

A justification of the procedure is in order. We first discuss the case 
in which at least one blocking constraint has been violated. To reduce 
the blocking probability for each queue whose constraint has been 
violated, a position must be added to this queue (by Property 4). The 
only other way to reduce the blocking probability is to add a server, 
but this branch has already been terminated. Assume that NS servers 
have been removed since the last ‘‘optimal” feasible solution and that 
the cost of all the positions added since the last ‘‘optimal’’ feasible 
solution is G. If G > NS X C, then the new parameters cannot give 
a lower cost solution and, by Property 4, no lower cost solution for 
this M exists. Therefore, we terminate the branch with M servers and 
begin the search of the branch with M — 1 servers (increment NS by 
1) with the present number of positions. At least this number of posi- 
tions must be considered, since, by Property 3, the reduction of M 
results in an increase in the blocking probabilities. If G < NS X C, 
this set of parameters may be a lower-cost solution. Therefore, we 
determine the system characteristics and see if the solution is feasible. 

We now prove that the procedure converges in a finite number of 
steps to a global optimum. 


Lemma 1: If a feasible solution for a given M has been found, the branch 
corresponding to that M can be terminated. 


Proof: Trivially, any further feasible solutions with that JZ must be 
more expensive, since these solutions must have more waiting positions. 
Q.E.D. 


Theorem 1: The algorithm terminates in a finite number of steps. 


Proof: First, we know that an initial feasible solution can be obtained 
in at most 


l + 
| ENP - Mo| 
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steps, as discussed earlier, where [2 ]* = max (0, x). Second, for any 
given M, the corresponding branch of the solution tree will be ter- 
minated in a finite number of steps. That is, since C is finite and C; is 
nonzero for all 7, then in a finite number of steps we will either find a 
feasible solution for this M, reach the point where the cost of the posi- 
tions added for this M exceeds C, or violate the delay constraint. In 
the first case, by Lemma 1, we know that we can terminate this 
branch. In the second case, since none of the positions can be removed 
and feasibility be maintained, no feasible solution of lower cost exists 
for this 7. In the latter case, the procedure terminates. The number of 
iterations performed for a given M is bounded by the number of times 
positions are added before the cost of these additions exceeds C. 
Finally, since the maximum number of servers that need be con- 
sidered is finite, we reach the case in which the delay constraint is 
violated in a finite number of steps (at most, M® values of M). 
Since there are only a finite number of values of M to be considered 
and since, for each M, only a finite number of steps are performed, the 
algorithm terminates in a finite number of iterations. Q.E.D. 


Theorem 2: The solution (N*, M*, Z*) obtained upon termination of the 
algorithm is a global optimum. 


Proof: Assume that another soneatatun (N, 7, Z) exists, such that 
all constraints of (I) are satisfied and Z < Z*. First, consider the case 
in which M@ > M*. By construction, the branch corresponding to 
must have been searched and, as indicated in Theorem 1, the branch 
would have been terminated in a finite number of steps. If this branch 
had produced a feasible solution with a lower cost, it would have been 
retained. Hence, this case is not possible. 

Second, consider the case in which M < M*. From the algorithm, 
we know that the termination of the procedure implies that the delay 
constraint has been violated. We therefore know that either the branch 
with M produced a feasible solution with a cost larger than Z* (or else 
it would have been retained), or M is smaller than the value of M 
when the procedure terminates. If the latter is true, then MZ cannot 
produce a feasible solution since the delay and blocking probability 
for (N*, M) must be greater than those for (N*, M*) by Property 3; 
and, to reduce this delay, positions would have to be removed that 
would result in at least one blocking constraint being violated. Con- 
sequently, this case is a contradiction and, hence, (N*, M*, Z*) is the 
optimum. Q.E.D. 

One point should be noted: For the higher levels of blocking (20.05), 
the solution (N, Mo) is generally a feasible solution. However, as 
a result of the reduction in offered load to the servers because of the 
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blocking, M (9) can be reduced before any constraints are violated. This 
is of importance since, in practice, such systems have been engineered 
in such a manner that the configuration is (N®, M(o)), which is ob- 
tained from use of the Erlang B and Erlang C formulas, as described . 
earlier. 

A more concise mathematical summary of the algorithm follows. 


Optimization Algorithm: 
Step 1: Initial Feasible Solution 


(t) Determine V{? (¢ = 1, 2, ---, 2) from the Erlang loss formula 
where NV{ = min {n|B(n, ai) S b;}. 

(it) Determine M,o) from the Erlang delay formula where 
Mo) = min {m|W(m, a) S d}. Let k = 1. 


kth iteration: 


(172) Calculate ByY(N©, M (-1), a) and D(N©), M (x-1), a). If the 
set of constraints (21) are satisfied, let M® = M x1. The 
initial “optimal” feasible solution is N* = N®, M* = M, 
and Z* = >°4_,CiN; + CM*. Set j = 1, 1N® = N%*, and go 
to Step 2. 


If at least one constraint of (21) is not satisfied, set Ma, 
= Mu-1) + 1, increment k, and return to (2772). 


Step 2: Solution Improvement 


In this step, the superscript 7 refers to the jth “optimal” value of 
M; for a given value of j, the subscript r refers to the rth value of N;. 


jth iteration: 


@) NS =1,G% =0,r = 1. 
(tt) Reduce number of servers by one, M% = MG) — 1. If 
M® = 0, go to Step 3. 
(dit) Calculate Bi(,NG™, M%, a) and D(,NG», M), a). 
(a) If DN‘, MM, a) > d, go to Step 3. 
(b) If (21) and (22) are satisfied and if G <= NS X C, then 
store (NN, M™) as the new “optimal” solution. That 
is, set N* = ,N&), M* = M® and 


T* = > C; NOY + OM, 
t=1 
Set iN‘ = ,N‘~, increment j, and go to (2). 
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(c) If (21) and (22) are satisfied but G“ > NS X C, then 
reduce the number of servers by subtracting 1 from M“), 
increment NS, and return to (i772). 

(d) If at least one blocking constraint of (21) is violated, i.e., 
SO = {7|Bi.NC, M®, a) > b;} is not empty, then add 
Dies C; to G™. Let 


riNf-) = >NVf-» a ¢ Sa 
and 
maNP- = NEP 41 1 8%, 


If G) > NS X C, then decrement M“) by 1, increment 
r, and go to (a2). If G < NS X C, increment r, and go 
to (272). 


Step 3: Termination 


If D(,.NG», M®, a) > d, then stop the procedure. The global 
optimum is (N*, M*, Z*). 


4.3 Numerical example 


To illustrate the algorithm, we examine a simple two-queue system 
that represents an automatic call distributor system used for credit 
checking. The company has subscribed to two Inward wats bands 
with a cost per trunk of $800 and $500. The two trunk groups each 
receive 15 erlangs of traffic in the busy hour. Calls, on the average, 
are 45 seconds in duration. The subscriber has requested a 5-second 
speed of answer and blocking objectives of P.10 and P.05, respectively. 
The monthly cost of an attendant is $750. With these parameters, 
we begin the algorithm by using the Erlang loss formula with 15 
erlangs and the delay formula with 30 erlangs to obtain (N™®, Mo), 
which are (18, 20; 34). As shown in Table I, this initial solution is feas- 
ible and hence will be stored as our tentative optimal solution, 
(N, MO), 

We proceed to Step 2 of the algorithm in an attempt to improve the 
initial feasible solution. By decreasing M® by one, MY = M® — 1, 
we obtain the system parameters (18, 20; 33) and the system char- 
acteristics (0.0912, 0.0504; 0.404), which indicate that this is not a 
feasible solution. Since the blocking constraint for trunk group 2 is 
violated, a trunk must be added to this group at a cost of $500. The 
accumulated cost of the additional trunks, $500, is less than the ac- 
cumulated cost of the removed attendants, $750. Therefore, we 
proceed by obtaining the system characteristics for this set of param- 
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Table |— An example of the optimization procedure 


Initial Feasible 


Solution 1] 18 20 34 | 0.0887 | 0.0481 | 0.201 | 49900 

2] 18 20 33 | 0.0912 | 0.0504 | 0.404 | 49150 

3 | 18 21 33 | 0.0926 | 0.0373 | 0.546 | 49650 

Sehation 4) 18 21 32 | 0.0970 | 0.0409 | 0.945 | 48900 
Tacievenient 5 | 18 21 31 | 0.1034 | 0.0460 | 1.550 | 48150 

Pp 6 | 19 21 30 | 0.0939 | 0.0573 | 2.990 | 48200 

*'7 | 19 22 30 | 0.0976 | 0.0466 | 3.422 | 48700 

8 | 19 22 29 | 0.1120 | 0.0572 | 4.976 | 47950 

Termination 9) 20 23 28 | 0.1220 | 0.0676 | 8.864 | 48500 


* Optimal solution. 
System Parameters: 


a, = 15 erlangs, bi = 0.10, ci = $800 
a2 = 15 erlangs, be = 0.05, co = $500 
HT = 45s, d=5s, C = $750. 


eters, i.e., (18, 21; 33). As Table I indicates, this is a feasible solution, 
is a cost improvement, and hence is stored as the tentative optimum. 

As shown in Table I, the procedure continues from this point until 
(20, 23; 28), at which point the delay constraint is violated. The 
optimal solution is (19, 22; 30) at a cost of $48,700. We should note 
that the parameters (20, 23; 29) were not examined since the accumu- 
lated cost of added trunks, $1300, was greater than the accumulated 
cost of removed attendants, $750. 

The above solution is the global optimum for this constrained 
problem. However, practitioners might suggest that (18, 21; 31) isa 
more realistic design since, in fact, the blocking constraint is “‘es- 
sentially” satisfied (0.1034 vs 0.1000). This can be incorporated in the 
design procedure by allowing any solution that is within ‘‘e’ of a 
blocking objective to be retained. The algorithm can then be applied 
as before. 


V. SUMMARY 


In this paper, an analysis of a particular multiserver, multiqueue 
service system has been presented. Examples of this type of system 
are the directory assistance systems used in the telephone companies 
and credit verification bureaus used by the credit-card industry, which 
use automatic call distributors. Expressions were derived for the 
equilibrium-state probabilities, and four properties of the system char- 
acteristics, overall average delay, and the blocking probabilities for each 
queue were given. 

These results were used in developing a procedure to obtain a least- 
cost system configuration to satisfy a given set of single-hour grade-of- 


622 THE BELL SYSTEM TECHNICAL JOURNAL, MARCH 1975 


service constraints. That is, the procedure determines the number of 
waiting positions for each queue and the number of servers required to 
satisfy constraints placed on blocking probabilities and average delay 
at minimum cost. The work reported here should form the basis for 
the development of a practical method of traffic engineering and 
administration for small automatic call distributor systems. 
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Blocking is considered for an N-trunk group of exponential servers 
with Poisson-offered load whose rate parameter varies with time. The 
infinite trunk case is solved by means of a rapidly convergent series of 
Poisson-Charlier polynomials. This solution is used to obtain practical 
approximations of blocking probability, transition probabilities, and 
recovery function for general time-variable offered load in the finite trunk- 
group case. An integral equation is derived satisfied by the blocking 
probability in the general case. In the situation of constant offered load, 
two additional methods are derived for providing easily computable 
approximations; one based on the integral equation, the other based on an 
approximate inversion of the Laplace transform. To aid in the latter 
approximation, bounds on the roots of Poisson-Charlier polynomials are 
obtained; in particular, an approximation ts obtained for the dominant 
root. The inversion of the integral equation is studied with the purpose of 
providing the basis for future investigations of errors of approximation. 
Curves are provided for a number of examples permitting comparison of 
exact and approximate solutions. 


I. INTRODUCTION 


The main purpose of this paper is to present a discussion of the 
time behavior of blocking in a fully available N-trunk group for any 
initial state with exponential servers and with Poisson-offered load 
whose rate parameter a(t) itself may be considered to vary with time; 
that is, the probability of 7 calls arriving in a time interval (0, ¢) is as- 
sumed given by 


: 

; I aude | 

exp (— f a(u)du) WA? +, 7.=-0, 1, 2,2. Cd) 
0 : 

The service rate is taken equal to 1, so that a(¢) is measured in erlangs. 

The problem of blocking with time-variable offered load was con- 

sidered by Palm! for finite trunk groups and by Khintchine? for infinite 
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trunk groups. The impetus for this is the need felt for more accurate 
computation of blocking probabilities and correlation information’ than 
can be obtained by quasi-stationary analyses, that is, by the use of 
equilibrium formulas in which the offered load parameter is replaced 
by its instantaneous value. Lack of statistical equilibrium renders 
this approach inaccurate. The time-variable aspect of the input stream 
should be carefully distinguished from other statistical descriptions 
such as peakedness,* since the effects on the system are separately 
identifiable. It has been reported, for example, that offered load and 
peakedness determinations from carried usage, peg count attempts 
offered to a group, and overflows are misled by the time variability 
of the offered load. 

Palm had proposed an interesting method of accounting for the time 
variability of the offered load, i.e., his “slow variations’ model. In 
this model, it was assumed that the actual ordering in time of a(é), 
that is, the functional dependence of a(t) on ¢ could be ignored if a(t) 
varied slowly. He replaced a(t) by a random variable with an incom- 
plete gamma function distribution. Thus, traffic functions such as 
blocking may be obtained from their equilibrium values by averaging 
over the appropriate gamma distribution. This model, however, re- 
quires further elaboration in view of the investigations of Iversen,’ 
who showed that the Palm approach does not correctly model the 
empirical data collected in the extensive Holbaek measurements of 
Danish telephone traffic. Iversen found that the correct time variation 
of the traffic could not be ignored. 

The trunk provisioning procedure whereby one uses the average 
offered load over a busy hour to achieve a required grade of service 
results, in some cases that were considered, in only a small under- 
estimate of the required number of trunks as calculated by the methods 
of this paper. Since the standard method is convenient, this may be 
viewed as substantiation of the approach. 

Essential for the methods of this paper is a Volterra integral equa- 
tion derived in Section II satisfied by the blocking probability, 
Py(t, N), experienced by a load of a(t) erlangs offered to an N-trunk 
group. Exact analytical solution of this equation is not useful, but 
numerical methods may be advantageously used. An important fea- 
ture of the equation, nonetheless, is that it permits studying errors of 
approximation and, in one instance (Appendix D), was directly used 
in the construction of an approximate operator for studying the 
transient response in the case of constant a(¢). Appendix A presents an 
explicit representation of Py(t, N) for general a(t) by means of an 
infinite Neumann expansion. Inequalities for Py(t, N) and truncation 
error estimates for the Neumann series are also given. 
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The infinite trunk group case forms the basis of the approximations 
developed in Section IV that are applicable to the general case of time 
variable a(¢). Although this case was solved by Khintchine,? Section III 
presents new representations in terms of rapidly convergent Poisson- 
Charlier expansions. Truncation error estimates are obtained, and the 
rate at which the state probabilities approach the Poisson form is 
assessed. To aid in the use of Poisson-Charlier polynomials, Appendix 
B provides a short discussion of their properties, especially providing 
convenient means of expanding a function into a Poisson-Charlier 
series. Since the state probabilities of the infinite trunk group system 
are often close to Poisson, this form of representation is very useful. 
The Poisson-Charlier expansion expresses the deviation of a function 
from the Poisson form. Further, in Section IV, the Poisson-Charlier 
polynomials are used to express the transition probabilities in explicit, 
closed form. 

The approximations of Section IV are applicable to time variable 
a(t), and are developed from the infinite trunk group solution by 
renormalization appropriate to the finite trunk group. To facilitate the 
use of the approximations, closed expressions are obtained for the 
relevant infinite trunk group solutions. This approximation procedure 
gives rise to the useful notion of a ‘‘modified offered load.’”’ One of the 
approximations obtained was, in fact, already obtained by Palm.! 
This approximation is particularly interesting because it uses the 
Erlang loss function, B(N, a), for which rapid methods of computation 
are avallable.*.® A special case of the approximations for transition 
probabilities is that for the recovery function, which is important in 
the discussion of correlation properties’ and, hence, in the determina- 
tion of variances of traffic parameter estimators. 

The constant offered load case is studied in Appendix C. Although 
the solution for the state probabilities is known,® the integral equation 
formulation appears to be new. Certain advantages are obtainable from 
this formulation. The errors of approximations to the state probabilities 
satisfy the same integral equation but with a different inhomogeneous 
term; thus, the more general integral equation is studied, leading to 
methods for investigating the quality of approximation. For this 
purpose, a natural Banach space (uniform norm over [0, ~]) is 
introduced, in which the integral operator is bounded and has a 
bounded inverse. Of course, the known Laplace transform of the 
transition probabilities is immediately obtained as a corollary. The 
integral equation is also used in Appendix D as a tool for the con- 
struction of an approximate solution (the scaling method) correspond- 
ing to an arbitrary initial state. Appendix C also presents several 
bounds on the required roots of Poisson-Charlier polynomials by 
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obtaining general bounds on the largest and smallest members of sets 
of positive numbers subject to isoperimetric constraints; the results 
are then specialized to the Poisson-Charlier polynomials. One of the 
bounds for the dominant root is explicit and easily calculable. Its 
accuracy appears to be good (see Fig. 8). 

Appendix D provides two approximations for the case of constant 
offered load. The scaling approximation, whose genesis was suggested 
by 8S. Horing, is constructed by obtaining an approximate time in- 
variant of the transition blocking probability from the initially empty 
state. The subsequent generalization of this approximation to arbi- 
trary initial states is then obtained by means of the integral equation. 
The Laplace transform approximation is constructed by an adaptation 
of the Widder formula’ for the inversion of the Laplace transform. It 
requires the determination of the dominant root; but, depending on 
the needs, it may be made arbitrarily accurate. It has the interesting 
property that, under certain conditions, it provides bounds for the 
exact solution. 

A discussion of some results and graphical illustrations is given in 
Section V. In testing the quality of the modified offered load approxi- 
mations, high change rates of a(¢) were chosen, in fact, much higher 
than would occur in practice. The errors of approximation increase with 
increasing rate of change of a(¢), hence, the examples chosen indicate 
much higher errors than one would expect to encounter in the practical 
application of these methods. 


Il. INTEGRAL EQUATION FOR BLOCKING 


It is the object of this section to establish Theorems 1 and 2, which 
provide integral representations of the binomial moments (13) of the 
probability distribution (2) of the number of busy trunks, and corollary 
1 of Theorem 2, in which an integral equation is given for the prob- 
ability that all trunks are busy at a given time. 

Let & = &(t, N) be the number of trunks busy at time ¢ in an N-trunk 
group, and P; = P;(t, N) the corresponding probability, 


PLEG, N) = j] = P(t, N). (2) 
The probability generating function g(t, ¢, N) is given by 
N 
g = Hetty) = ¥" Pst, NDS. (3) 
j=0 
At the point ¢ + dt, one has 
gi + di,§,N) =9 + ag, E(t -+ dt, N) = &€+ dé (4) 
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and, hence, 


g + dg = Ey {gtEl¢4| Eq}, (5) 


in which L[¢*| £] is the conditional expectation of ¢¢ given £, and FE; 
is the expectation over the probability distribution of £ The boundary 
at 7 = N necessitates a further analysis of (5). One has 


g+ dg = (1 — Py)E, {Stee | < NJ} 
+ Py, {gES*|— = N]}. (6) 
The probability distribution of dé is 
Pldé = —l1jé=j,0S jf <N]= jdt, 
Pldt=0/F = 37,05 7 < NJ =1- (a+ jd, 
P[dé= 1 = j,0 Sj < N] = ad, (7) 
P[dt = 0|—E = N] =1—Ndt, 
Pldé = —1|/& = N] = Ndi; 
hence, 
ELS#|é= j7,0S 7 < NJ =14 (af —a— jt je“)dt, 
E[s#|—= NJ =1—N(1— 5“)dt. 
From (6) and (8), one has 


(8) 


“4 = (1— Py)E,[¢t(at —a—é+ YE < NY 
— PyE,[&Oo( —- Dl|E=N], (9) 
a = (1 — Py)a(é — DEE E< NJ - ¢& — DEALER); 
hence, 
+e 1) i = a(t — Ig — a(t — DE*Pw. (10) 


The infinite trunk group does not require the analysis of (6) nor the 
boundary conditions ( = N) of (7), hence (9) becomes 


0 
ap = as — Dg — & — DEE AI (11) 
Thus, the corresponding equation for the infinite trunk group is 
og ny eee 
+ GD oe = alt — Do. (12) 
The binomial moments 8,(é, NV) are defined by 


g.(, NV) = & Px) (2), (18) 
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in which (2) designates the binomial function defined by 


(jd eee, 


s s! 


(J) =1. 


The binomial moment generating function l(t, w, NV) = disdo Bs(t, N)w* 
is given by 


IV 


s2 1, 


(14) 


L(t, w, N) = gt, 1+ w, N). (15) 
Thus, the differential equation satisfied by I is 
al abe . 
ah + wa awl — aw(1 + w)*Py. (16) 
The corresponding equation for the infinite trunk group is 
al al 
ai +w oo awl. (17) 


Equation (16) is a linear partial differential equation that can be solved 
by the following device (method of characteristics). Let @ be a new, 
independent variable and set 


l= 1(6), 
w= w(8), (18) 
t = ¢(6). 


Then, comparison of 
dl aldt al dw 


do aidd * ow do a” 
with (16) yields the set of equations 
ee aw — aw(1 + w)*Py 
dé ae 
dw 
dt 
de 


whose solution for J is then obtained. To exhibit the solution conveni- 
ently, let 


Bea) =e | * era(s)ds, (21) 


630 THE BELL SYSTEM TECHNICAL JOURNAL, MARCH 1975 


and 
A = A(t) = a(t, 0) = e-*xa(t), (22) 


where * designates convolution product over (0, ¢). The solution of 
(16) takes the form 


L(t, w, N) = 1(O, we-#, N)eAv 
— wf ert] + wet*]Na(r)Pyw(t, N)dr. (238) 
0 


Similarly, the solution of (17) is 
L(t, w, ©) = 1(0, we’, ~)edr. (24) 


To obtain the binomial moments themselves, it is convenient to intro- 
duce the Volterra operator K., 


t 
Kf = [ Kfar, (25) 
0 
defined by the kernel 
sal qi N 
= pall —(s—J) (t-1) 
K.(t, 2) Balgsi as) DeMa(r). (26) 
Since the Laguerre polynomial L& (—<)! is given by 
@/_») . < a n+ta 
Le(-2) = Ea"), (27) 


the kernel K,(t, 7) may be written more compactly; thus, 
K,(t, 7) = a(rye #9 LAST) (—aet-*), (28) 


The following theorems may now be stated. 


Theorem 1: The binomial moments, B;(t, ©), for the infinite trunk group 
are given by 
8 . Ad 
B3(t, 00 ) = 2 B_;(0, eG yt 
7=0 J 
Proof: The coefficient of w* in the expansion of the right-hand side 
of (24) yields the result. 


Theorem 2: The binomial moments, 8,(t, N), for the finite trunk group 
satisfy 
B. (6, N) = Bs(t, 0 ) = K.Py. 


Proof: The coefficient of w* in the expansion of the right-hand side of 
(23) provides the required result. 
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Corollary 1: The probability, Pn(t, N), that all trunks are busy satisfies 
the integral equation 


Py(t, N) = By(t, ~) — KyPy, 
ain which 


“a ,, Ad 
By (i, 00 ) = a By_;(0, Wye tre 
j=0 | 


Ky(t, 7) = a(rye“X 9 L® 1 (—aet-’). 
Proof: For the finite trunk group (13) shows that 
By(t, N) = Pyit, N). 


Hence, the integral equation follows from Theorem 2. The explicit 
expressions for By(t, ©) and Ky(é, r) are obtained from Theorem 1 and 
(28), respectively. 

The special case of all trunks free initially leads to a somewhat 
simpler integral equation for Py(t, N). This is given in the following 
corollary. 


Corollary 2: When all trunks are oi free, Pn(t, N) satisfies 
Py(t, N) = Ni7 KyPy. 


Proof: The initial probability distribution P;(0, N) has the form 
P;0,N)=1 7=9, 


29 
= 0 j > 0. (29) 
Hence, the binomial moments satisfy 
(0, N) =1 = 0, 
6.(0, N) 3 ey 


= 0 s> 0. 
The equation for By(t, ©) given in the first corollary now yields 
AN 
The result of the corollary follows. 
The probability Py(t, N) corresponding to all trunks initially busy 


is called the recovery function; it satisfies the following integral 
equation. 


Corollary 3: When all trunks are initially busy, one has 


Py, N) = e~Nt*Ln(—e'A) — KyPy. 
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Proof: We have 


P;(0, N) = 0 Ox<j<N 
At ) ) f = Jj < } (32) 
= 1 Z=N. 
Hence, the required binomial moments are 
N 
8.(0, V) = (7) ): (33) 


The result follows from the equation for By(t, ©) and from (27). 
A noteworthy case occurs when a is constant, then the integral 
equation for Py(t, N) becomes of convolution type. Thus, let 


Ky(t) = ae-™*L®_,[—a(et — 1)]. (34) 
Then one has Corollary 4. 
Corollary 4: When the offered load a is constant, Py(t, N) satisfies 
Py(t, N) = Bv(t, <0 ) — Ky*Py. 
Proof: Use of (21) and K y(t, 7) as given in Corollary 1 yields the result. 
For constant a, an equilibrium distribution P;(«, NV) exists.? Let 
x(a) = 2 (35) 
Then 


P,(, N) = O5j8N. (36) 


a’ 
j!Sy(a)’ 
The notation B(N, a) is used for the blocking probability Py(«, N) 
and is referred to as Erlang’s loss formula.® Corresponding to the 
equilibrium distribution, one has the binomial moments 6,(0, N) and, 


hence, the moments for the infinite trunk group given in Theorem 1. 
These moments will be denoted by 6{(t, a). It may be noted that 


lim ps(¢, a) = &- (37) 
to 0 §! 


The integral of Ky(t) that is useful in error analyses may now be 
easily obtained. 


Theorem 3: When the offered load is constant, we have 


; _ BY, a) 
i Ku(a)dr = Beye gy —h 


I Keds = Sey 4, 





in which Ky(r) is gwen in (84). 
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Proof: Since B(N, a) is the equilibrium solution of the integral equa- 
tion of Corollary 4, and since 6%(t, a) corresponds to the equilibrium 
state, the solution of the integral equation is constant and equal to 
B(N, a); thus, 

B(N, a) + Ky*B(N, a) = Balt, a). 


This immediately implies the first equation of the theorem. The second 
follows on considering t >, and using (36) and (37). 

It may be useful to observe that the positive character of the general 
operator Ky in Corollary 1 immediately supplies the inequalities 


By(t, ©) — KwBw(t, ©) < P(t, N) < Bui, @). (38) 


A Neumann-series solution of the integral equation of Corollary 1 is 
discussed in Appendix A. Higher-order inequalities of type (38) are 
also given. 


Ill. THE INFINITE TRUNK SOLUTION 


It will be convenient to express the solution for P,(t, ©) in terms of 
Poisson-Charlier polynomials"! whose relevant properties are dis- 
cussed in Appendix B. The probability distribution of the number of 
busy trunks for the infinite trunk case was considered by Khintchine? 
and, for constant offered load, by Karlin and McGregor.’ Theorem 4 
presents a rapidly convergent form of the solution in terms of Poisson- 
Charlier polynomials valid for any initial state. This solution will 
be the main tool for the construction of approximations to distribu- 
tions in the finite trunk case. 

From (15), let 1(0, w, ©) be given by 


1(0, w, ©) = % 8:(0, co wi, (39) 


then the binomial moment generating function for the infinite trunk 
case is, from (24), 


I(t, w, ©) = eb ¥ Bi(0, -o)e-wi, (40) 


The mean, pn, of this distribution is the coefficient of w; hence, 
B= A at Boe, (41) 


in which po is the mean of the initial distribution. One may now state 
Theorem 4. 
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Theorem 4: The probability distribution, P(t, ©), of the number of busy 
trunks in an infinite trunk group has the convergent representation 


P,(t, ©) 


I 


Y(a, L) s e-*d;G,(a, bt), | 





J —] v 
dj; Dm = u’B;-»(0, 0), 


v=0 


provided that 
1(0, 2 — 1, ©) = 4 B,(0, »)(2 — Ii = EPO, ~)z! 
j=0 j=0 
converges for |z| <r(r > 2). 


Proof: Use of (40) in (101) which the choice \ = yp as given in (41) 
yields 








Cw) = emo" B,(0, co ew, (42) 
7=0 
Hence, 
— pit : (1) v 
oj = ei YS ys, (0, «) (43) 
and, from (95), 
Pall, 2) = W(t, #) & eGi(e, w). (44) 
= 
Let 
oe OED ng. 
d; zs es t HoB;-»(0, 0), (45) 
»=0 ov! 
then 
cj = etd, (46) 


and the formula of the theorem follows. A theorem of Uspensky' 
states that the general representation of (95) is valid in the sense that 
the series converges to F(x) if the series }>7-» F(x)z* has radius of 
convergence greater than 2. Since 


L(t, 2 — 1, oo ) ia Eat Az 0, (z => le“, 0 |], (47) 


the radius of convergence of [(¢, z — 1, ©) is greater than 2 by the 
condition on 1(0,z2— 1, ©) stated in the theorem. Hence, by 
Uspensky’s theorem, the representation of (44) is valid for all ¢ 2 0. 

A truncation error estimate for the series representation of Theorem 
4 is given in Theorem 5. 


NONSTATIONARY BLOCKING 635 


Theorem 6: 


|P.(t, ©) — ¥(2, u) — ¥(e, n) > e-#d,G;(, u) | 


2 2 k+1 e—(kt+l)t 
= Sd rp nw 
rm (=) 1 — (2/w)e-* CO Ue) 


k=1, R>w>0O, t>inZ, 
in which R ts the radius of convergence of 1(0, w, ©). 


Proof: Since 


lw(a, B)G; (a, b) | = |p (a, B) | < 24, j 2 0, (48) 
we have 
a e~*ld;| |b(x, u)G;(x, w)| S zu e~ #23 d;|. (49) 


Also, from (45), 
|dj| = 3 Bj-+(0, 2) S erorl(0, w, oo )w~ (50) 
Thus, 
& e~*ld;| |W (a, u)G;(a, uw) | S exml(0, w, ©) & e~2iw-F (51) 


k+1 —(k+1)t 
=) € (52) 


< ghow = Se a ay 
SOTHO, thy 2) (- 1 — (2/w)e 
The conditions of the theorem ensure the convergence of I(0, w, ~) 
and of the series of (51). 

The corollaries below follow directly under the conditions of Theorem 
5. 


Corollary 1: 
k 
Pz(t, ©) = W(x, w) De e-M*djG;(a, w) + O(e-*Y 4), 
j=0 
Corollary 2: 
P < 4 a now] (0 
|P:(t, ») — (a, w)| S wl Gwe? (0, w, ©). 


Corollary 3: 
Pz(t, ©) = w(x, w) + O(e**). 


Thus, the distribution quickly becomes nearly Poisson with the time 
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variable parameter yu, regardless of the initial distribution. In fact, if the 
initial distribution is Poisson, one might anticipate P,(t, ©) to remain 
Poisson for all t = 0. This is asserted by Theorem 6. 


Theorem 6: If P,(0, ©) is Potsson with parameter yo, then 
P,(t, ©) = (a, u). 


Proof: The binomial moments, 8,(0, ©), are 


B.0,0) =, 820. (53) 
It follows, from (45), that 


wh (i ; 
a= (7)\(-yy=0, j>0. (54) 
J ev»=0 \ PV 
The result is now obtained from Theorem 3. 


IV. APPROXIMATIONS 


The Neumann series (75) is an explicit solution of the integral 
equation of blocking given in Corollary 1, Theorem 2. For constant 
offered load a, the resolvent kernel solution (108) and Theorem 11 are 
also available; however, especially when WN is large, these solutions 
are not convenient. It is therefore important to have approximations 
that lend themselves to computation for large N in a sufficiently 
convenient manner. Three such approximations have been developed, 
namely: the ‘‘modified offered load’ approximation that is useful and 
fairly accurate in the general case, that is, for time variable offered 
load, the “scaling” approximation, and the ‘Laplace transform” 
approximation, which are applicable only to transient phenomena 
under constant offered load. The scaling approximation is also fairly 
accurate and does not require factorization of polynomials. The Laplace 
transform approximation consists in fact of an infinite set of approxi- 
mations of arbitrarily high accuracy. It usually requires finding a 
single root—the so-called dominant root—of an appropriate poly- 
nomial. The scaling and Laplace transform approximations are dis- 
cussed in Appendix D. Appendix C provides approximations for the 
required dominant root. The modified offered load approximation 
is presented below. 

Let P;,,(t, N) be the probability that the N-trunk group started 
from state 7 at time 0 and proceeded to state x at time #, then P;,.(¢, ©) 
may be computed from Theorem 4 using 


8,0, «) = (4): (55) 
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An approximation, P;,,(t, N), for P;,.(t, N) is given by 
Py s(t, 00 ) . 


Pi2(t, N) == 
y. P(t, co ) 


(56) 


This approximation is suggested by the following considerations. One 
has 
lim P;,.(t, N) = Pi,2(t, ©). (57) 
N->o 


Hence, the approximation should be accurate even for time-varying 
offered load when N is large. Furthermore, when a is constant, since 
limi. = a, one has 


sone a 
lim Piz(t, N) = Seay? (58) 


which, as previously indicated, is the exact equilibrium distribution; 
hence, the approximation should be accurate when ¢ is large even for 
time-varying a(t), provided d(é) is small. Since, by the law of total 
probability, 


Pall, N) = ¥ Ps(0, N)Piclb, N), (59) 


one can construct an approximation, P,(t, N), to P.(t, N) by use of 
P;,.(t, N); thus, 


P,(t, N) = ¥ P.(0, N)P;,2(t, N). (60) 


To facilitate the use of (60), Theorem 7 expresses P,,,(¢, ©) in finite 
form. 


Theorem 7: 
P;.2(t, ©) = (1 — e~*)teA "ati, —(et — 1)A]. 


Proof: From Theorem 4 and (55), we have 





Pi,2(t, ©) = ¥(2, w) ps e-Hd,G;(2, n); (61) 
i (=1)’. 
a= yo? Pie) (62) 


Comparison of (62) with (90) shows that 


WP ee Pet 
d; = qi! G;(4, 1). (63) 
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Let 








c= = * Gui, a, (64) 
then 
OGG Ati, Bee 
C(w) = Y -—* GG, 2)w, (65) 
j=o 7! 
which, by comparison with (89), may be rewritten 
C(w) = e-*e-™(1 + e-*tw)*. (66) 
Let 
B(w) = e#°C(w) = e4¥(1 + e-'w)* (67) 
and , 
g(z) = B(z — 1) = e~4e42(1 — e-* 4 e-tz)*. (68) 
Then 
P3,2(t, 00 ) 


= ene E (FEUER DAL 6g 


=0\%— V v 


Hence, comparison of (69) with (90) finally yields 





Pi,2(t, ©) = (1 — e-*jie-4 Mo, —(e'—1)A]. (70) 


Immediate corollaries are the following: 


Corollary 1: 
Pall, @) = eA 9" ¥ po, @)(1 — GL, —(e — DAD 


Corollary 2: 


P;,.(t, N) = G.[i, —(et — 1)A] 


> (x !/v!)(—A)’?*GL[1, — (et — 1)A] 


Of particular interest are the functions Po,n(t, N) and Py,w(t, N); 
the first describes the progression of the system from initially empty to 
blocked, and the second describes the recovery of the system from an 
initially blocked condition to the blocked condition again. The latter 
function is called the “recovery function.”’’? The following formulas 
are obtained from Corollary 2. 


Pow(t, N) = B(N, u). (71) 


The general principle of approximation employed, namely, the re- 
normalization of an appropriate solution for the infinite trunk group 
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case, allows one to state, by use of Theorem 6, that, whenever one starts 
from an approximately Poisson initial state, an approximation to 
Py(t, N) is 

Py(t, N) = BN, p). (72) 


This approximation was known to C. Palm.! The parameter y is re- 
garded as a modified offered load. 
The recovery function approximation obtained from Corollary 2 is 


Gw[N, —(et — 1A] 


Py w(t,N) =7y 
> (N !/v!)(—A)’""GLLN, — (et — 1A] 


(73) 


V. NUMERICAL EXAMPLES 


Tor the purpose of providing some idea of the accuracy of the ap- 
proximations developed in Section IV and Appendix D, curves were 
drawn up comparing exact and approximate solutions for a group of 
ten trunks. These curves illustrate nonstationary behavior. Figure 1 
shows the scaling and modified offered load approximations for a step 
input problem in which a = 7 erlangs is offered to an initially empty 
group. The scaling approximation of (166) was used, and (71) was used 
for the modified offered load approximation. Apparently, for this situa- 
tion, the scaling approximation is somewhat more accurate. 


0.10 
0.09 
0.08 
0.07 
z= 0.06 
= 
co 
< 0.05 
je) 
a 
= 0.04 TRUNK GROUP = 10 TRUNKS 
OFFERED LOAD a=7 ERLANGS 
0.03 B (10, 7) = 0.07874 (EQUILIBRIUM LEVEL) 
——— _ EXACT PROBABILITY 
0.02 — — —— SCALING APPROXIMATION 


—--—— MODIFIED OFFERED LOAD APPROXIMATION 





0 1 2 3 4 5 6 7 
t (TIME) 


Fig. 1—All trunks empty initially—scaling and modified offered load approximations. 
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TRUNK GROUP = 10 TRUNKS 
OFFERED LOAD a= 7 ERLANGS 
B (10, 7) = 0.07874 (EQUILIBRIUM LEVEL) 
EXACT PROBABILITY 
— —~ ——— LAPLACE TRANSFORM APPROXIMATION 


wae em MODIFIED OFFERED LOAD APPROXIMATION 


P (PROBABILITY) 





t (TIME) 


Fig. 2—Recovery function—Laplace transform and modified offered load approx- 
imations. 


The recovery function approximations of (73) and (187) are com- 
pared to the exact solution of (185) in Fig. 2. The approximations are 
correct at the extremes ¢ = 0,¢ = o, and track the exact curve reason- 
ably well. The approximation of (73) is more accurate initially and is, 
of course, also applicable when the offered load is time-variable; how- 


TRUNK GROUP = 10 TRUNKS 
OFFERED LOAD a=4+ 4t— 1.6t2,0<t<5 
=4,5<t<10 


TRUNKS IN EQUILIBRIUM INITIALLY UNDER 
OFFERED LOAD a= 4 ERLANGS 


EXACT PROBABILITY 


———— MODIFIED OFFERED LOAD 
APPROXIMATION 


P (PROBABILITY) 





t (TIME) 
Fig. 3—Pulse response—modified offered load approximation. 
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TRUNK GROUP = 10 TRUNKS 
OFFERED LOAD a=7+3S!N 47t ERLANGS 
EXACT PROBABILITY 

—— — — MODIFIED OFFERED LOAD APPROXIMATION 


P (PROBABILITY) 





t (TIME) 


Fig. 4—All trunks empty initially—modified offered load approximation. 


0.08 


TRUNK GROUP = 10 TRUNKS 
OFFERED LOAD a= 5+2SIN 0.057rt 
EXACT PROBABILITY 
—————— MODIFIED OFFERED LOAD APPROXIMATION 


0.07 


P (PROBABILITY) 


0.06 


0.05 





5 6 7 8 
t (TIME) 


Fig. 5—All trunks empty initially—modified offered load approximation. 
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P (PROBABILITY) 


P (PROBABILITY) 


0.07 


0.05 


0.04 


0.03 


0.06 


0.05 


0.04 


0.02 


TRUNK GROUP = 10 TRUNKS 
OFFERED LOAD a=5+2SINO.17t 
EXACT PROBABILITY 
——-—— MODIFIED OFFERED LOAD APPROXIMATION 


5 6 7 8 
t (TIME) 


Fig. 6—All trunks empty initially—modified offered load approximation. 


TRUNK GROUP = 10 TRUNKS 
OFFERED LOAD a=5+2SINO0.57t 
EXACT PROBABILITY 


———-— MODIFIED OFFERED LOAD 
APPROXIMATION 


t (TIME) 


Fig. 7—All trunks empty initially—modified offered load approximation. 
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ever, the approximation of (187) is simply the case n = 0 of Theorem 
16. Considerable enhancement of accuracy may, for example, be ob- 
tained by using n = 1 or even higher values of n. 

The modified offered load approximation, (72), is compared to a 
pulse response in Fig. 3. It is seen that, despite the rapid variation of 
a(t) (as high as 14 erlangs/call duration) and the large range of prob- 

‘ability values, the approximation well imitates the course of the 
response. 

Figures 4, 5, 6, and 7 illustrate the modified offered load approxi- 
mation applied to sinusoidal inputs. In all cases, the trunk group was 
initially empty. Figure 4 shows the response, starting from ¢ = 0, to 
a=7+3sin4at. This may be considered to have wide excursions 
compared to the constant term 7, and rapid oscillations, i.e., the period, 
T, is 5. The exact curve is seen to be well imitated by the 
approximation. 

One may consider the total error to consist of two components, an 
evanescent part arising from the specific initial state and a component 


TRUNK GROUP = 8 TRUNKS 
EXACT ROOT 
—— — — APPROXIMATION 





0 1 2 3 4 5 6 7 8 
a — OFFERED LOAD 


Fig. 8—Upper bound approximation to dominant root. 
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resulting from the rate of change of the offered load itself, that is, a 
function of d(¢) which vanishes when a(t) = 0. Figures 5, 6, and 7 
are intended to illustrate the latter component above; hence, the time 
scale starts at five. The periods of 40, 20, and 4, respectively, were 
chosen to reflect the effect of d(¢) on the approximations. To provide 
clearer comparison, the probability scales of the graphs have been 
expanded. 

The lower bound of Theorem 13 provides an upper bound on the 
dominant root. A comparison with the exact values for an eight-trunk 
group taken from Bene?’ is given in Fig. 8. 


VI. NEEDED INVESTIGATIONS 


Much work remains to be done to provide a satisfying and viable tool 
for fully available trunk group analyses. One may mention error estima- 
tion of the approximations suggested in Section IV and Appendix D, 
the investigation of new approximations, such as studying the con- 
sequences of using a refined scaling approximation or an improved 
modified offered load, and the study of (J + Ky) in the Banach 
space X, for variable a(t). This, in turn, would permit new approxi- 
mations to be constructed and would provide improved means of error 
investigation. Relatively little is known about the behavior of the 
zeros of Poisson-Charlier polynomials, especially in the present 
context, as functions of a, N. It is hoped this paper will provide an 
impetus for further investigation. 
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APPENDIX A 
Neumann-Series Solution 
The Neumann-series solution of the integral equation for blocking, 
Py(t, N) = By(t, ©) — KwPw, (74) 
is 
Py(t, N) = Bult, ©) — KyBw(t, ©) + KiB (t, ©) — +--+, (75) 


which, of course, is convergent for all ¢. The positivity of Ky implies 
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the system of inequalities 


2k+1 ; 2k : 
d (—1)'KyBw(t, ©) < Py(t, N) < 2 (—1)’Kxbr(t, ©), k 2 0, 
j=0 j=0 


(76) 
which generalize (38). 
Truncation error estimates for (75) will now be obtained. The in- 
equalities (76) yield 


—K¥By(t, ©) < Pw(t, N) — ¥ (-I5Kw, ©) <0, (77) 


0 < Py(t, N) — 3 (—1)iKiBy(t, ©) < K#By(t, 0). (78) 
7=0 


Hence, it is only necessary to bound K%t'8y(t, ©). Let 


a(t) SG, (79) 

then, from (27) and (84), 
Ky(t, 7; a) S Ky(t— 17; 4@) S Ae, (80) 
A = GL$_\(—@). (81) 


The dependence on a is explicitly shown in (80). One similarly obtains 
2 N Gi 
By(t, ae a) $B= 2 Bw—s(0, N) j! ) t2 0. (82) 
3 = . 


One may now state Theorem 8. 


Theorem 8: 
At 
nBu(t, ©) S pio oe 
Proof: One has 
Ky(t; a) S Ae. (83) 
Hence, 
r—1 
Ky,,-(t; G) < At (r— 1)! e', (84) 


in which Ky,,(t; @) is the r-fold convolution of Ky(t; @) with itself. 
Convolving this with By(t, ©) finally yields 


4" 


Ky,+(t; @)*By(t, ©) S$ B—>- (85) 
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APPENDIX B 
Poisson-Charlier Polynomials 


Some properties of Poisson-Charlier polynomials" are developed in 
this appendix, especially with a view to convenient representation of 
functions by series expansion. Let 


W(n,) =X, pe AO xa, (86) 


Then the polynomials G;(x, 4) are defined by 


di 


dvs ¥(x, ») = W(z, A) = Gi(a, \)W(a, A). (87) 


The Taylor expansion 
W@rA+D = E Tyo,» (88) 
yields the generating function 
et (1 + +) FG (89) 
rN 7=0 gj! 
Thus, explicit formulas for G;(x, \) are 


G;(z, \) = ty (-1( = )s 


j-—vjyv! 


o = (—)P (7) yiar(2). (90) 


The first few polynomials are 


Go(a, ) => 1, 
Gila, ») =; (@—), 
1 (91) 
G2(x, A) = rc [a2 ~ (24 + l)x + d?7], 
G3(x, A) = < Ca! — 3(A + 1a? + (8d? + 3d + 2)a — AF]. 
A recurrence relation derived from (89) is 
Gix(a, ) = ==2=*Gy(,.) -264@r). 92) 
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The following inner product is defined for functions f(x), g(x) of the 
discrete variable x: 


(f, 9) = Le, YI @)9(@). (93) 


The Poisson-Charlier polynomials are orthogonal with respect to this 
inner product 


(G;, Gr) = 0, j zk 


(94) 
(Gi, @) = 


Accordingly, the coefficients c; in the expansion of a function f(z) 
in the form 


(2) = W(@, ) X esGj(2, ») (95) 
are given by 
oj = FE ile, NIC). (96) 


For the purpose of the present investigation, a more convenient mode 
of determining c; is achieved by obtaining their generating function; 
that is, 


Cw) = ¥ eww! (97) 
From (96), one has 
Cw) = % sa) XPS? Gila, ». (98) 
Hence, from (89), 
Cw) = ¥ f(e)(1 + w)*. (99) 
Let 
B(w) = ¥ fe) + w)* (100) 


Then B(w) is the binomial moment generating function of f(x) and 
C(w) = e>"B(w). (101) 


From (101), the first few coefficients c; are obtained in terms of the 
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corresponding binomial moments 8; of f(x); thus, 
Co = Bo, 
C1 = B1 — dBo, 
Cz = Bz — ABi + ZA*Bo, 
C3 = Bs — AB2 + Zd7B1 — GA*Bo. 
A useful choice of the parameter \ is suggested by (102), namely, 


(102) 


A= B1/Bo, (103) 
which yields c; = 0 and 
2 
Co = Bo — a ’ 
BiB Bi = 
2 1 
Cs = Bs — “a 1 age 


For a probability distribution, one has 68> = 1, and the choice (103) 
for \ implies that \ is equated to the mean of the distribution. In this 
case, one has 

C2 = 3(0? — 4), 

C3 = g(a — 3o7 + 2Qy), 


in which yp is the mean, o? the variance, and e@ the third moment about 
the mean of the distribution. 


(105) 


APPENDIX C 
Constant Offered Load—Dominant Roots 
In this appendix, the integral equation for the constant a case, 


namely, 
Py(t, N) = Brit, ©) — Ky*Pwit, N), 


Ky(t) = ae-"'LY_[—a(et — 1)], (106) 
Bw(t, ©) = eM! = By_;(0, N) ee 


is studied. In fact, the somewhat more general hn 


f@® + Kw*fl) = 9), (107) 


is resolved. This presents a considerable advantage over the solution 
of (106), since the errors of approximations to Py(t, N) satisfy (107), 
and hence may be studied by means of Theorem 10 and its corollary. 
Solutions for blocking have been obtained in the literature,’ but do 
not provide a means for error analysis. For the practical utilization 
of the solutions, bounds for the exponents occurring in the explicit 
representation of the resolvent kernel will also be obtained. 
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Integral equation theory” asserts that a resolvent kernel, Qy(t), 
exists with the property 


f® = 9g) + Qn«g(). (108) 
Laplace transformation will be used to study (107) and (108). 


Theorem 9: The Laplace transforms, Kw(s), On(s), of Kn(t), Qn(Q), 
respectively, are 


(—1)%a%Gw(—s — 1, a) 


Kxn(s) = (ao ee 
Cy(s) = GB ++(8 + N) = (=1)"a"Gu(=s — 1, a), 
al (=1)*a"G@n(—s — I, a) 
Proof: One has, from (27), 
L®_(- ae 2 (aa) a (109) 
and hence, ae 
Ky(t) = ae-™! x ¢ : 1 i (ef — 1)%. (110) 
Thus, 


ReGrse ps ty i )e [re e-Wte-I(L — etidt, (111) 


Letting X = e~‘, one has 


Ry) = aE (7, 


The integral in (112) is the beta function, B(N +s — j, 7+ 1). 
Hence, 


a f} ; 
) A i XNA — X)ide. (112) 


One has the following transformations: 
Rulo= % (5) weaweesoaD cu) 
w= yea & (4 eet y--G+N-A, (15) 
Rul) = ey ee (GUD GHD, — 16) 


RvO= Gay eem BOG) ee( FG) an 
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Thus, the required formula for Ky(s) is obtained from (117) by 
comparison with (90). From (107), one has 








ft Ky+f =@ (118) 
Hence, 4 
a Z ade sehen Je: 
f= tak, 7 27-74% (119) 
thus, : 
hey An Ky(s) ; 
Qn(s) i+k© (120) 


This, together with the formula for Kn(s), yields the required expres- 
sion for Qy(s). 


Corollary: 
G:( —&, a) 


P;,n(s, N) = a ee Sa), 


Proof: Use of Theorem 2, Corollary 4, Theorem 9, and eq. (178) with 
An expression equivalent to this corollary was given by Takacs. 
It will be useful now to introduce the Banach space, X, of functions 

f(t) that are bounded and measurable over (0, ~) and normed by 


fll = sup [f@|. (121) 
t20 
One may now state Theorem 10. 


Theorem 10: The operators I+ Ky and (I+ Ky)™ are bounded; 
further, 
|Z + Kw = Sy(a), 


N+ Kyl = 1+ f° l@w@lae 


an which the operator norms are those induced by (121). 


Proof: The quantity ||ZJ + Ky]|| is obtained directly from Theorem 3 
and the formula 


|Z + Kwll = 1+ f° Kw(@at. (122) 
Since the polynomials Gy(x, a) are orthogonal over (0, ~), it follows 
that the zeros, Pi, ---, Py are distinct, real, and positive; hence, the 


Zeros, 71, °++, Tv, Of Gy(—s — 1, a) are distinct, real, and less than 
minus one. The Paley-Wiener theorem? applied to Oy(s) now asserts 
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that J + Ky has a bounded inverse with norm C given by 


C= 0+ Kyl =14 [° lQv@lat (123) 
An immediate corollary is the following 


Corollary: The integral equation 


ft+Ky*f= g, 


gE X, 
possesses a solution f € X satisfying 


fll S Cllgtl, 
fim | f(t}| < C lim |g(| 
t-»00 td 
In particular, if mise g(t) = 0, then limi... f(t) = O 
Proof: The result follows from 


f= 9+ Qn*g 
and the definition of C. 


(124) 
The following theorem provides a representation of Qy(t) and an 
estimate of C. 


Theorem 11: 


N II (rj eh v) 
Q(t) = > 
II 





1 " 
erit C < » ——s 
I (y—1) [rol UL bry — rol 
Ser: ar, 
Proof: One has 
(—1)%a%Gy(—s — 1, a) = (8 — 11)--- (8 — ry) (125) 
and, hence, from Theorem (9), 
A; (sd) eee PNY — (sr) Ses — Fy) 
= Ss 126 
Qn(s) (s — 11)--- (8 — Tw) ( ) 
The partial fraction expression for Oy(s) is 
N il (ry v) 1 
On(s) = % —_ __. (127) 
ia II (7; mz ry) 
1 


we 
Kl 
& 
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Thus, 





N 
N II (rj v) 
Qy(t) = —) erit, (128) 
j=l I (1; — t;) 
yay 
Also, one has - 
N IL [rj + >| 
l\QvO| S Dw erit (129) 
7 ~ih [rj = ry | 
aa 
and, accordingly, : 
: IT |r +» 
C=1+f" lQvoldts1+ a ey - (130) 
. 7dr TT [ry — rl 


rr 


To make the results of Theorem 11 accessible for estimation, partic- 
ularly in numerical applications, and for the Laplace transform ap- 
proximation of Appendix D, bounds for the roots, 7;, will be obtained. 
The generating function, g(z), for the equilibrium probability distribu- 
tion, P;(«, N), of (86) may be written as 


_ Sw(az) 





g(z) = Svar (131) 

Thus, the mean, m, and variance, o?, are 
m = ali — B(N, a)], (132) 
o = m— (N — mjaB(N, a). (1338) 


We now have Theorem 12. 
Theorem 12: (Benes) 11 > — (m/o?). 

To obtain further bounds, the following lemmas are needed. 
Lemma 1: py 2 ++: 2 pi > O, 


pitts tpn =S1,  pit-:: +py= 82, pitspn = D, 


then 





IA 


< Sit VN — 1)(Ns2 — 83) 
N ) 


pyv& 


IA 
IIA 


p Pl 


in which p ts the small positive root of 
e(s1 — p)* 1 = D(W — 1)*7. 


NONSTATIONARY BLOCKING 653 


These bounds are sharp. 


Proof: The equations 


pit-:: tpv=S, (134) 
pit -++ + py = 8 (135) 
may be written 
pit -+: + py-1 = 81 — pn, (136) 
pi t+ +++ + py—1 = S82 — py. (137) 
Let 
py = Max py (138) 


over all allowable sequences. Then the sum (1387) is minimum when py 
is replaced by py. This occurs, however, only when 


2 * 


Pl = *°° = pn-i = wo . (139) 





Thus, from (137), 





me BND 
(W—1) (SoH) = 5 — of (140) 
The solution of (140) for py is 


si + V(N — 1)(Nsz — 8%) 
ee (141) 


* 
Pu = 
This proves the upper bound of the lemma. The inequality is attained 
for the vector (p1, ---, pw) defined by (139) and (141). 
Similarly, one may write 


pe + -++ + pw = 81 — pi, (142) 
p2-* pn = D/p. (1438) 

Let 
pi = min py (144) 


over all allowable sequences. Then the product (143) is maximum when 
px is replaced by p;. This occurs only when 





yee ap 


Thus, 





Sfp et rn 
| met | = Dist (146) 
and the equation of the lemma defining p has been established. The 
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inequality is attained for the vector (1, ---, pv) defined by (145) and 
(146). 


Lemma 2: py 2 ++: 2 pi >0O 


Then 


in which p is the small positive root of 
1 N-1 
D (s — *) = p(N — 1)", 


The bound ts sharp. 


Proof: One may write 


1 1 1 
—+----+—=81--, (147) 
ps PN pi 
p2-* pn = D/px. (148) 
Let ” : 


over all allowable sequences. Then the product (148) is maximum 
when p; is replaced by p;. This occurs when 


N-1 





p, = +++ = py = i (150) 
S_1 — 
Py 
Thus, 
N-—1)*> 
1 | = D/p; (151) 
S_1 — 7 
Py 





and the equation of the lemma is established. The inequality is at- 
tained for the vector (p1, ---, pw) defined by (150) and (151). 

The application of the lemmas to the polynomials Gy(z, a) is ac- 
complished by identifying p1, ---, pw with its zeros. For this purpose, 
the form of Gy(az, a) given in (90) will be recast, by the help of the 
Stirling numbers of first kind," into standard form; that is, 


Gy(a, a) = > an nt™. (152) 


m= 
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The Stirling numbers, S7, are defined by 


,—1 . ; 
Tl (@ —») = ¥ Ste". (153) 
v=0 m=1 
Thus, one has 
N N—m N 
Gy(e, a) = (-)* + Yam ZS (—1rar-w (1) sp, (154) 
The sums s:, s2 are accordingly given by 


sy = (2) +N, (155) 


n= st—6(7)- @a+4) (3) -2#(7): (156) 


The reciprocal polynomial, that is, the polynomial whose zeros are 
pi’, ***, pw’ is given by 


N 
X*Gy(a7, a)Gw(a, a) = YS Gne™. (157) 
m=0 
Hence, the analogous quantities s_; and s_, = pp? + --- + py? are 
given by 
a | 
= > ge ares (158) 
v=] 
5 xd ma | 
s2=s,—-2 > -NMa’ oD = (159) 
y=2 V j=l ae 
in which 


N® = 1, N®™ = N(N —1)---(N—v+ JD), y>0O. (160) 
The upper bound of Lemma 1 now establishes Theorem 13. 


Theorem 18: 


N 

see ate ne er ee i 

mt VN = Nese) oo 
-at+VW— Wa — a) 
=~ —_-W 





Proof: The upper bound is immediate. The lower bound is obtained 
by applying the upper bound of Lemma | to the reciprocal equation. 

A numerical illustration of Theorem 13 is provided by the zeros of 
Gio(x, 7) used in obtaining the recovery function plotted in Fig. 2. 
They are 0.332811, 2.05847, 4.06653, 6.31227, 8.81308, 11.5197, 
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14.6407, 18.0255, 22.0872, and 27.1438. The lower and upper bounds 
given by the theorem are 0.32964 and 36.82292, respectively. It appears 
that the lower bound may well be usable as an approximation to 1. 
The accuracy of this when used to approximate —ri = 1 + pi is 
illustrated in Fig. 8. The exact values of —7ri are taken from Benes.’ 
This provides an upper bound for 71 which, together with the available 
lower bounds, is useful for error investigations. 


APPENDIX D 


Approximations—Constant Offered Load 


It was suggested by 8. Horing that an appropriate scaling between 
Py(t, N) and P,(t, k) may exist; that is, a function F[Pw(t, N), 
P,(t, k)_] may exist, which would be approximately independent of ¢ 
and which may, therefore, permit the approximate determination of 
Py(t, N) in terms of P(t, 1), for example, thus permitting large trunk 
groups to be studied in terms of the behavior of small ones. Since it is 
feasible to use (108) for small trunk groups, this would constitute an 
approximation of the solution of (106) for large trunk groups. 

Consider the following 


B(kl, a) _ Lk! Si(a/k)* 
BG, a/b) = GT Sula) (161) 


The ratio S:(a/k)*/S,.(a) is approximately independent of a since 
Si(a/k)* = Sii(a) & e* for k large. Thus, 


B(kl, a) 


BU, a/kyt oe) 


is approximately independent of a. It would seem, therefore, that the 
ratio (162) is approximately a time invariant of (106), especially for 
large t; that is, the function 


Priilt, kl; a) 


Pilt, [o/h ee 
is approximately equal to the ratio (162) when ¢ is large; thus, 
ae Pit, 1; a/k) lf 
Pyil(t, 1; a) & B(kl, a) | Beary a/b) (164) 
We have, from (107) and Theorem 11, 
: = a — pe-(atl)t 
Pit, 1; a) rere {l-—e } (165) 


when the trunk is initially empty. Hence, from (164) with / = 1, 
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k = N, 
Py(t, N) 2 BIN, a){1 — er Welrnieyy (166) 


for the case when all N trunks are initially empty. 
It is desired now to extend the approximate solution (166) for 
application to any initial condition. Define qn (t) by 


qn (é) = B(N, a) {1 — eI nyy, (167) 
Then a convolution operator with kernel Ly (t) will be constructed so 
that one has 


AN 
qn + Ly*Qn = = 


i? (168) 


exactly. The operator Ly will then be taken to approximate the opera- 
tor Ky of (106). Accordingly, an approximation Py to Py correspond- 
ing to initial states other than all trunks empty is defined by the 
equation 

Py + Ly#Py = By. (169) 


Theorem 14: The Laplace transform, Ly(s), of the approximating 
Ly (t) 1s 


ror (+ @ayeit}) 


rete + Dr omyzi) 


a s + »L(a/N) + 1] 
-(#+1)" Sy (a) yy Stee) ET). 


eae (4 +1) Sr) 


Proof: One has 
dv = B(N, a) is em] — e (aly) Wat, (170) 
0 


Let u = [(a/N) + 1]t, then 


.~ _  B(N, a) 
iv (G/N) +1 


The substitution x = e~“ yields 


B(N, a) 


e~ [si (aI N+1)]u(] — e~“) Ndu. (171) 
0 


ee Gi i, * ple! (alN-H)1-1(] — a) dz; (172) 
thus, 
s 
iy = 2O® wl oar] (173) 
LS ae 1| 
(a/N) +1 
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One also has 


eo N!T(s) 
—st — p—-t\)\N = “ 
i e(L — Md = Be eT (174) 
The Laplace transform of each term of (168) now yields 
2 i) ae 


Use of (173) yields the result of the theorem. 
It may be noted that Ly has an impulsive component at ¢ = 0 whose 
value is 


i —N 

(F +4 1) Slayer: (176) 
For large N, this is nearly zero. The effect of this is to create an 
error at ¢ = 0; that is, Py (0, N) ¥ Py(0, N) if Bv(0, N) + 0. The 


larger N is (for fixed a), the smaller the discrepancy. 
Equation (169) is studied in Theorem 15. 


Theorem 18: The solution of (169) is 


5 Nt . | GIN) = al 
(a/N) +1 

¥ 8i(0, N)a-T(s + 4), 
5 _ _ TD +3) _d 
Py = p> B;(0, N)a “F(D). qn(t), D= di 


Proof: From Corollary 1, Theorem 2, one has 


N. ai 
By = X Bx-s(0, Nae ae say ir (177) 
=> . 
hence, ¥ ( ’ 
. _, Fist j 
= N J), 
Transforming the terms of (169) yields 
Py(1 + Ly) = By, (179) 


and hence the result of the theorem is obtained from Theorem 14 and 
eq. (178). The inversion of the transform, Py, by use of differentiation 
follows on use of (173) in Py. 
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Another method of approximation useful for constant offered load is 
based on an approximate inversion of the Laplace transform. Let 


7G)= i “e-ef(uldu, s>0; (180) 
then : 
LOD™ got Fio (9 = =f e~*u" f(u)du. (181) 


The function (s*t!/n!)e-*“u" is a probability density function for 
any s > 0, n 2 0 whose mean is (n + 1)/s and variance (n + 1)/s?. 
Letting s = (n + 1)/t, the mean and variance are ¢ and #/(n + 1), 
respectively; hence, Korovkin’s theorem on sequences of positive 
functionals!® establishes the Widder inversion formula’: 


eS) ca ae. 
lim i. +1 F(™) (3) 


nao . 


= f@. (182) 





s=(n+1)/t 


The above discussion forms the basis for Theorem 16. 


Theorem 16: For ¢ > O, let the transform of e*f(t), namely, f(s — ©), 
exist for s > O, and let e«'f(t) be convex in t > O, then 


f) st I ponfon(s — 6 n20, t>0. 





s=(n+1)/t : 
Proof: Jensen’s inequality applied to (181) in the form 


es srt fim(s — 6) = i i esures” f(u)du (183) 
n\ n! Jo 

establishes the theorem. By virtue of (182), when similarly modified 

for the function e'f(t), the dexter of the inequality always provides 


an approximation to f(t) even when e*'f(t) is not convex. 
Corollary: f(s — 6) exists for s > 0, and e*f(t) is convex in t > 0 
implies 
1 3/1 
fi) Spey - ep t> 0. 
Proof: The case n = 0 of Theorem 16. 


If f(s) should have a dominant pole, it is usually advantageous to 
choose e equal to the negative of that pole. 

The above corollary will now be applied to obtaining an inequality 
for the recovery function. The corollary to Theorem 9 shows that 


Gv(—s,@) 
sGy(—s — 1, a) 





Py.v(s, N) => (184) 
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The inversion of Py,yv(s, N) is readily accomplished; the result is’ 





Py.w(t, N) = BOW, a) - & =  (1-—>—). cass) 
j=l 15 ij ri Ts 
It follows from (185) that 
e-"' Py v(t, N) — B(N, a) ] (186) 


is convex for ¢ > 0 and that the corresponding Laplace transform exists 
for s > 0, hence, by the above corollary, one obtains 


Py, v(t, N) s BIN, a) 


4 gt Gy[-(1/t) — n, a] B(N, a) 


Tt tal ee | 087 
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